guide.tex

\documentclass[titlestyle=hang,11pt]{elegantbook}

\author{Jinxiong \& Zhang}
\email{jinxiongzhang@qq.com}

%\zhtitle{优美的\LaTeX{} 书籍}
%\zhend{模板}
\entitle{A Guide to Data Science}
%\enend{Template}
\version{3.00}
\myquote{From mathematics to algorithms.}
\logo{logo.png}
\cover{cover.pdf}

%green color
   \definecolor{main1}{RGB}{0,120,2}
   \definecolor{second1}{RGB}{230,90,7}
   \definecolor{third1}{RGB}{0,160,152}
%cyan color
   \definecolor{main2}{RGB}{0,175,152}
   \definecolor{second2}{RGB}{239,126,30}
   \definecolor{third2}{RGB}{120,8,13}
%blue color
   \definecolor{main3}{RGB}{20,50,104}
   \definecolor{second3}{RGB}{180,50,131}
   \definecolor{third3}{RGB}{7,127,128}

\usepackage{makecell}
\usepackage{lipsum}
\usepackage{texnames}
\usepackage{xcolor}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{graphicx}
\usepackage[colorlinks,linkcolor=blue]{hyperref}
\usepackage{amsmath}
\usepackage{animate}
\usepackage{subfigure}

\begin{document}
\maketitle
\tableofcontents
\mainmatter

\chapter{Abstract}

This is a brief introduction to data science. \href{https://www.wikiwand.com/en/Data_science}{Data science} can be seen as an extension of statistics.

\begin{itemize}
    \item $\textcolor{red}{\text{ All sciences are, in the abstract, mathematics.}}$
	\item $\textcolor{red}{\text{ All judgments are, in their rationale, statistics.  by C.P. Rao}}$
\end{itemize}

\href{http://norvig.com/chomsky.html}{On Chomsky and the Two Cultures of Statistical Learning} or \href{https://projecteuclid.org/euclid.ss/1009213726}{Statistical Modeling: The Two Cultures (with comments and a rejoinder by the author)} provides more discussion of statistical model or data science.

\includegraphics[width = .8\textwidth]{john-von-neumanns-quotes-7.jpg}

\chapter{Statistical Foundation}

"Statistics, in a nutshell, is a discipline that studies the best ways of dealing with randomness, or more precisely and broadly, variation", \href{https://statistics.fas.harvard.edu/about}{Professor Xiao-Li Meng said}.

\begin{table}[h]
	\centering
	\caption{On mathematical model}
	\begin{tabular}{|c|}
		\hline
		\includegraphics[width = .95\textwidth]{Albert.jpg} \\
		\includegraphics[width = .95\textwidth]{model.jpg} \\
		%\hline
		%\includegraphics[width = .8\textwidth]{john-von-neumanns-quotes-7.jpg}
	\end{tabular}
\end{table}




\section{Basic Probability}
Probability theory is regarded as the theoretical foundation of statistics,which provides the ideal models to measure
the stochastic phenomenon, randomness, chance or luck.
It is the distribution theory bridging probability and statistics.

The project \href{https://seeing-theory.brown.edu/#secondPage/chapter5}{Seeing theory} is a brief introduction to statistics and probability theory
and provides a \href{https://seeing-theory.brown.edu/doc/seeing-theory.pdf}{.pdf file}.


\subsection{Random Variable}

\begin{definition}{Random Variable}{r.v.}
 A random variable is a mapping
 \[ \color{green}{X:\Omega \rightarrow \mathbb{R}} \]
that assigns a real number $X(\omega)$ to each outcome $\omega$.
\begin{itemize}
 \item $\Omega$ is the sample space.
 \item $\omega$ in $\Omega$ are called sample outcomes or realization.
 \item Subsets of $\Omega$ are called $\mathrm{Event}$.
\end{itemize}
\end{definition}

\begin{definition}{Probability}{Pr}
A function $P$ is called probability function with respect to the known sample space $\Omega$ and all the event $\forall \omega \in \Omega$ if
\begin{center}
\begin{tabular}{|l|c|}
 \hline
Properties& Conditions\\
 \hline
Nonnegative & {$P(\omega) \geq 0 \forall \omega \in \Omega$}   \\
 \hline
Normalized & $P(\Omega) = 1$  \\
 \hline
Additive & \text{If $\omega_1, \omega_2, \cdots \in \Omega$ are disjoint, then $P(\bigcup_{i=1}^{\infty}\omega_i)=\sum_{i=1}^{\infty}P(\omega_i)$}  \\
 \hline
\end{tabular}
\end{center}

\end{definition}

These properties are called \href{https://en.wikipedia.org/wiki/Probability_axioms}{Kolmogorov axioms}.

\subsection{Discrete Distribution}

Random variable or event is countable or listed, then we can list the probability of every event.
If so, we can compute the probability of a random variable less than the given value. For example, $P(X\leq 10)=\sum_{i=1}^{10}P(X=i)$ if the random variable $X$ only takes positive integers.
\begin{definition}{Probability mass function}{pmf}
The probability mass function  of a discrete random variable $X$ is defined as
$$f_X(x)=P_X(X=x)$$
for all $x$.
\end{definition}
The left subindex $\mathrm{X}$ is to point out that the probability is  with respect to the random variable $X$ and it can omit if no confusion.

\subsection{Continuous Distribution}

If random variables may take on a continuous range of values, it is hard or valueless to compute the probability of a given value. It is usually to find the cumulative distribution function.

\begin{definition}{Cumulative Distribution Function}{CDF}
 A function $F(x)$ is called cumulative distribution function (CDF in short) if
\begin{itemize}
 \item $\lim_{x\to -\infty}F(x) = 0$ and $\lim_{x\to \infty}F(x) = 1$ ;
 \item The function $F(x)$ is monotone increasing function with $x$;
 \item $\lim_{x\to x_0^+ } F(x) = F(x_0)$ for all $x_0$.
\end{itemize}
\end{definition}
In probability, $\forall x$, $F(x) = P(X\leq x)$ for some random variable $X$.

For example, if the random variable $X$ has the cumulative distribution function (CDF)
\[ F_X(x)=\begin{cases} x, & \text{$x \in [0, 1]$}\\1, &\text{$x\in (1, \infty)$}\\0, & \text{otherwise}\end{cases}, \]
we say the random variable $X$ is uniformly distributed in $[0, 1]$.
We can see \href{https://www.encyclopediaofmath.org/index.php/Continuous_distribution}{this link}.


The probability density function can be seen as the counterpart of pmf.
For the continuous random variable $X$, its CDF $F_X(x)$ is not required to be differentiable.
But if so, we can compute the derivative of the CDF $F_X(x)$: $\frac{\mathrm{d}F_X(x)}{\mathrm{d}x} = f_X(x)$.

\begin{definition}{Probability Density Function}{pdf}
We call the function $f_X(x)$ is the probability density function with respect to the continuous random variable $X$ if
\[ F_X(x) = \int_{-\infty}^{x}f_X(t)\mathrm{d}t \]
for all $x$, where $F_X(x)$ is the probability cumulative function of the random variable $X$.
\end{definition}

For example, the random variable $X$ uniformly distributed in $[0, 1]$ has the probability density function
\[ f_X(x)=\begin{cases} 1, & \text{$x \in [0, 1]$}\\0, & \text{otherwise}\end{cases}. \]

If the pdf $f_X$ is positive in the set $\mathbf{S}$, we call the set $\mathbf{S}$ is the support set or support of the distribution.

One common probability - $\color{teal}{\text{normal or Gaussian distribution}}$ - is often given in the pdf:
\[ f(x|\sigma^2, \mu)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}, \]
where $\sigma^2$ is given positive and finite, $\mu$ is given and finite. And its support is $(-\infty,\infty)$.
\begin{center}
\includegraphics[width = 16cm]{Gaussian.jpg}
\end{center}

The density function of a power law distribution is in the form of
\[ f_X(x)= x^{-\alpha}, x\in [1,\infty) \]
where  $\alpha$ is real. It is also called \href{https://www.wikiwand.com/en/Pareto_distribution}{Pareto-type distribution}.

The power law can be used to describe a phenomenon where a small number of items is clustered at the top of a distribution (or at the bottom), taking up $95\%$ of the resources. In other words, it implies a small amount of occurrences is common, while larger occurrences are rare.

Every probability distribution is an ideal mathematical model that describes some real event.
Here is [a field guide to continuous distribution](http://threeplusone.com/FieldGuide.pdf).

\subsection{Bivariate Distribution}

The cumulative probability function and probability density function are really functions with some required properties. We wonder the bivariate functions in probability theory.
\begin{definition}{Joint Distribution Function}{j.d.f.}
The joint distribution function $F:\mathbb{R}^{2}\rightarrow [0,1]$, where $X$ and $Y$ are discrete variables, is given by
\[ F_{(X,Y)}(x,y)=P(X\leq x, Y\leq y). \]
\end{definition}
Their joint mass function  $f:\mathbb{R}^{2}\rightarrow [0,1]$ is given by
\[ f_{(X,Y)}(x,y)=P(X=x, Y=y). \]

\begin{definition}{Joint Density  Function}{J.D.F}
The joint distribution function $F_{(X,Y)}:\mathbb{R}^{2}\rightarrow [0,1]$, where $X$ and > $Y$ are continuous variables, is given by $F_{(X,Y)}(x,y)=P(X\leq x, Y\leq y)$. And their joint density > function if  $f_{X,Y}(x, y)$ satisfies that
\[ F_{(X,Y)}(x,y)=\int_{v=-\infty}^{y}\int_{u=-\infty}^{x}f_{(X,Y)}(u,v)\mathrm{d}u\mathrm{d}v \]
for each $x,y \in \mathbb{R}$.
\end{definition}


\begin{definition}{Marginal Distribution}{Marginal Dis}
The marginal distributed of $X$ and $Y$ are
\[ F_X(x)=P(X\leq x)=F_{(X,Y)}(x,\infty)\quad\text{and}\quad F_Y(y)=P(X\leq x)=F_{(X,Y)}  (\infty, y) \]
for discrete random variables;
\[ F_X(x)=\int_{-\infty}^{x}(\int_{\mathbb{R}}f_{(X,Y)}(u, y)\mathrm{d}y)\mathrm{d}u \quad\text{and}\quad  F_Y(y)=\int_{-\infty}^{y}(\int_{\mathbb{R}}f_{(X,Y)}(x,v)\mathrm{d}v)\mathrm{d}x \]
for continuous random variables and the marginal probability density function is
\[ f_X(x)=\int_{\mathbb{R}}f_{(X,Y)}(x, y)\mathrm{d}y,\quad f_Y(y)=\int_{\mathbb{R}}f_{(X,Y)}(x, y)\mathrm{d}x.\]
\end{definition}

Two random variables $\color{red}{X}$ and $\color{blue}{Y}$ are called as identically distributed if $P(x\in A)=P(y\in A)$ for any $A \subset \mathbb{R}$.

\section{Representation of Random Variable}

The random variable is nearly no sense without its probability distribution. We can choose the CDF or pdf to depict the probability, which depends on the case.

We know that a function of random variable is also a random variable.
For example, supposing that $\color{teal}{X}$ is a random variable, the function $\color{teal}{g(X)= e^{X}}$ is a random variable as well as $\color{teal}{g(X)=\ln(X)}$.

$\color{teal}{\text{The function of random variable must have its distribution}}$.
We can take it to generate more distributions.
\begin{theorem}
Suppose that $X$  is with the cumulative distribution function $\color{red}{\text{CDF}}$ $F_X(x)$ and $Y=F_X(X)$, then $Y$ is $\color{teal}{uniformly}$ distributed in the interval $[0, 1]$.
\end{theorem}
Let $X$ and $Y$ be two different random variables, $Z= X+Y$ or $Z = X\cdot Y$ are typical functions of random variables, especially $X$ is discrete while $Y$ is continuous.
%More on the Wikipedia page [random variable](https://www.wikiwand.com/en/Random_variable)

\subsection{Mixture Representation}

Let $P_1,P_2,\cdots, P_n$ be probability distribution and $\sum_{i=1}^{n}\lambda_i=1,\lambda_i > 0 \quad \forall i\in\{1,2,\cdots,n\}$, we can get
\[ P=\sum_{i=1}^{n}\lambda_iP_i \]
is also a probability distribution.
If $X$ is  distributed in $P$, i.e.$X\sim P$, its probability is computed as
\[ P(X=x)=\sum_{i=1}^{n}\lambda_iP_i(X=x) \]
for discrete random variable or
\[ P(X\leq x)=\sum_{i=1}^{n}\lambda_iP_i(X\leq x) \]
for continuous random variable.

 A random variable $X$ is said to have a mixture distribution if the distribution of $X$ depends on a quantity that also has a distribution.


Sometimes, it is quite difficult if we directly generate a random vector $X\sim f_X(x)$, but the augmented vector $(X, Z)
^⊤\sim f_{(X,Z)}(x, z)$ is relatively easy to generate such as [Box-Muller algorithm].
And we can represent $f_X(x)$ as the marginal probability distribution of $f_{(X,Y)}(x,y)$ in integral form
\[ \int_{\mathbb{R}}f_{(X,Z)}(x,z) \mathrm{d}z. \]


Thanks to Professor Tian Guoliang(Gary) in SUSTech, who brought me to computational statistics.
More information on probability can be founded in \href{http://probweb.berkeley.edu/probweb.html}{The Probability web}.

\section{Bayes' Theorem}

Bayes theorem is the foundation of Bayesian statistics in honor of the statistician \href{https://en.wikipedia.org/wiki/Thomas_Bayes}{Thomas Bayes}.

\subsection{Conditional Probability}

If the set $A$ is the subset of $B$, i.e. $A\subset B$, we know that if $x\in A$, $x\in B$.
We can say if $A$ happened, the event $B$ must have happened. However, if $B$ happened, what is probability when $A$ happened? It should be larger than the probability of $A$ when $B$ did not happened and it is also larger than the probability of $A$ when the events including $B$ happened. We know that $A\bigcap B=A$ when $A$ is the subset of $B$.
\begin{definition}{Conditional Probability}{Conditional Prob}
	Supposing $A, B\in \Omega$ and $P(B) > 0$, the conditional probability of $A$ given $B$(denoted as $P(A|B))$ is
	\[ P(A|B)=\frac{P(A\bigcap B)}{P(B)}. \]
\end{definition}
The vertical bar $|$ means that the right one is given or known and is parameter rather than variable.

Here is \href{http://setosa.io/conditional/}{a visual explanation}.

\begin{definition}{Statistical Independence}{Independence}
	If $P(A\bigcap B)= P(A)P(B)$, we call the event $A$ and $B$ are statistically independent.
\end{definition}

From the definition of conditional probability, it is obvious that:
\begin{itemize}
	\item $P(B)=\frac{P(A\bigcap B)}{P(A|B)}$;
	\item $P(A\bigcap B)=P(A|B)P(B)$;
	\item $P(A\bigcap B)=P(B|A)P(A)$.
\end{itemize}
Thus
\begin{itemize}
	\item $P(A|B)P(B)=P(B|A)P(A)$;
	\item $P(A|B)=\frac{P(B|A)P(A)}{P(B)}=P(B|A) \frac{P(A)}{P(B)}$;
	\item $P(B)=\frac{P(B|A)P(A)}{P(A|B}=\frac{P(B|A)}{P(A|B)}P(A)$.
\end{itemize}

\begin{definition}{Conditional Distribution Function}{Condintional Distribution Function}
	The conditional distribution function of $Y$ given $X=x$ is the function $F_{Y|X}(y|x)$ given by
	 \[ F_{Y|X}(y|x)=\int_{-\infty}^{x}\frac{f_{(X,Y)}(x,v)}{f_X(x)}\mathrm{d}v=\frac{\int_{-\infty}^{x}f_{(X,Y)}(x,v)\mathrm{d}v}{f_{X}(x)} \]
	 for the support of $f_X(x)$ and the conditional probability density function of $F_{Y|X}$, written $f_{Y|X}(y|x)$, is given by
	 \[ f_{Y|X}(y|x)=\frac{f_{(X,Y)}(x,y)}{f_X(x)} \]
	 for any $x$ such that $f_X(x)>0$.
\end{definition}

\begin{definition}{Conditional Probability for Continuous Random Variable}{Conditional Probability for Continuous Random Variable}
	If $X$ and $Y$ are non-generate and jointly continuous random variables with density $f_{X,Y}(x, y)$ then,  if $B$ has positive measure,
	\[ P(X\in A|Y \in B)=\frac{\int_{y\in B}\int_{x\in \color{red}{A}}f_{X,Y}(x,y)\mathbb{d}x\mathbb{d}y}{\int_{y\in B}\int_{x\in \color{red}{\mathbb{R}}}f_{X,Y}(x,y)\mathbb{d}x\mathbb{d}y}. \]
\end{definition}

\begin{definition}{Chain Rule of Probability}{Chain Rule of Probability}
\[ P(x^{(1)},\cdots, x^{(n)})=P(x^{(1)})\prod_{i=2}^{n}P(x^{i}|x^{(1),\cdots, x^{(i-1)}}) \]
\end{definition}

\begin{definition}{Conditional independence}{Conditional independence}
The event $A$ and $B$ are \href{https://www.wikiwand.com/en/Conditional_independence}{conditionally independent} given $C$ if and only if
\[ P(A\bigcap B|C)=P(A|C)P(B|C).\]
\end{definition}

\begin{theorem}{Total Probability Theorem}{Total Probability Theorem}
Supposing the events $A_1, A_2, \cdots$ are disjoint in the sample space $\Omega$ and $\bigcup_{i=1}^{\infty}A_i = \Omega$, $B$ is any subset of $\Omega$, we have
\[ P(B)=\sum_{i=1}^{\infty}P(B\bigcap A_i)=\sum_{i=1}^{\infty}P(B|A_i)P(A_i). \]
\end{theorem}

\begin{note}
\begin{itemize}
 \item For the discrete random variable, the total probability theorem tells that any event can be decomposed to the basic event in the event space.
 \item For the continuous random variable, this  theorems means that $$\int_{\mathbb{B}}f_X(x)\mathrm{d}x=\sum_{i=1}^{\infty}\int_{\mathbb{B}\bigcap\mathbb{A_i}}f_X(x)\mathbb{d}x.$$
\end{itemize}
\end{note}

\subsection{Bayes's Formula and Inverse Bayes' Formula}

\begin{theorem}{Bayes's Formula}{Bayes's Formula}
Supposing the events $A_1, A_2, \cdots$ are disjoint in the sample space $\Omega$ and $\bigcup_{i=1}^{\infty}A_i = \Omega$, $B$ is any subset of $\Omega$, we have
\[ P(A_i|B)=\frac{P(B|A_i)P(A_i)}{\sum_{i=1}^{\infty}P(B|A_i)P(A_i)}=\frac{P(B|A_i)P(A_i)}{P(B)}. \]
\end{theorem}

\begin{note}
\begin{itemize}
 \item The probability $P(A_i)$ is called prior probability of event $A_i$ and the conditional probability $P(A_i|B)$ is called posterior probability of the event $A_i$ in Bayesian statistics.
 \item For any $A_i$, $P(A_i|B)\propto\, P(B|A_i)P(A_i)$ and $P(B)$ is the normalization constant as the sum of $P(B|A_i)P(A_i)$.
\end{itemize}
\end{note}
\begin{center}
\includegraphics[width = .8\textwidth]{thomas_bayes.jpg}
\end{center}

The joint pdf of the random variables $X$ and $Y$ can be expressed as
\[ f_{(x,y)}(x,y)=f_{X|Y}(x|y)f_{Y}(y)=f_{Y|X}(y|x)f_{X}(x). \]
in some proper conditions.
Thus we get by division
$$ f_{Y}(y)=\frac{f_{Y|X}(y|x)f_{X}(x)}{f_{X|Y}(x|y)}. \eqno{(1)} $$

Integrating this identity with respect to $y$ on support of $f_{Y}(y)$, we immediately
have the \textbf{point-wise formula} as shown below
$$ f_{X}(x)=\{\int_{f_{X|Y}(x|y)\neq {0}}\frac{f_{Y|X}(y|x)}{f_{X|Y}(x|y)} \mathrm{d}y\}^{-1}.\eqno{(2)} $$
Now substitute (2) into (1), we obtain the dual form of IBF for $f_Y(y)$ and hence by symmetry we obtain the \textbf{function-wise formula} of $f_Y(y)$ at $y_0$ as shown in (f1), or the
sampling formula in (f2) when the normalizing constant is omitted.
\begin{align}
f_{X}(x)
&=\{\int_{f_{Y|X}(y|X)\neq {0}}\frac{f_{X|Y}(x|y_0)}{f_{Y|X}(y_0|x)}\mathrm{d}x\}^{-1}\frac{f_{X|Y}(x|y_0)}{f_{Y|X}(y_0|x)} \tag{f1} \\
&\propto \frac{f_{X|Y}(x|y_0)}{f_{Y|X}(y_0|x)}  \tag{f2}
\end{align}

\begin{table}[h]
	\centering
	\caption{The Inventor of IBF}
	\begin{tabular}{|c|}
		\hline
		\includegraphics[width = .75\textwidth]{kaing.jpg} \\
	\end{tabular}
\end{table}

There are more information on \textbf{inverse Bayes' formula}:
\begin{itemize}
		\item \url{http://web.hku.hk/~kaing/Background.pdf}
		\item \url{http://101.96.10.64/web.hku.hk/~kaing/Section1_3.pdf}
        \item \url{http://web.hku.hk/~kaing/HKSSinterview.pdf}
		\item \url{http://web.hku.hk/~kaing/HKSStalk.pdf}
\end{itemize}


\section{What determines a distribution?}

The cumulative density function (CDF) or probability density function(pdf) is roughly equal to random variable, i.e. one random variable is always attached with CDF or pdf. However, what we observed is not the variable itself but its realization or sample.
\begin{itemize}
  \item In theory, what properties do some CDFs or pdfs hold? For example, are they all integrable?
  \item In practice, the basic question is how to determine the CDF or pdf if we only observed some samples?
\end{itemize}

\subsection{Expectation, Variance and Entropy}

Expectation and variance are two important factors that characterize the random variable.

The expectation of a discrete random variable is defined as the weighted sum.
\begin{definition}{Mean}{Mean}
For a discrete random variable $X$, of which the range is ${x_1,\cdots, x_{\infty}}$, its mean $\mathbb{E}(X)$ is defined as  $\mathbb{E}(X)=\sum_{i=1}^{\infty}x_i P(X=x_i)$.
\end{definition}

The mean of a discrete random variable is a weighted sum. Average is one special expectation with equiprob, i.e. $P(X=x_1)=P(X=x_2) = \cdots = P(X = x_i) = \cdots \forall i$.

The expectation of a continuous random variable is defined as one special integration.
\begin{definition}{Expectation}{E}
 For a continuous random variable $X$ with the pdf $f(x)$, if the integration $\int_{-\infty}^{\infty}|x|f(x)\mathrm{d}x$ exists, the value $\int_{-\infty}^{\infty} xf(x)\mathrm{d}x$  is called the (mathematical) \textbf{expectation} of $X$, which is often denoted as $\mathbb{E}(X)$.
\end{definition}

Note that $\color{lime}{NOT}$ all random variables have expectation. For example, the expectation of standard \href{https://www.wikiwand.com/en/Cauchy_distribution}{Cauchy distribution}
$$f_X(x)=\frac{1}{\pi(1+x^2)}$$
is undefined.

The expectation is the center of the probability density function in many cases. What is more,
$$\arg\max_{b}\mathbb{E}(X-b)^2=\mathbb{E}(X)$$
which implies that $\mathbb{E}(X)$ is the most likeliest to appear in expectation.

\begin{definition}{Expectation of Random vector}{Exp}
 Now suppose that X is a random vector $X = (X_1, \cdots, X_d )$, where $X_j$, $j = 1, \cdots, d,$ is real-valued. Then $\mathbb{E}(X)$ is simply the vector
 $$(\mathbb{E}(X_1), \cdots , \mathbb{E}(X_d))$$
where $\mathbb{E}(X_i)=\int_{\mathbb{R}}x_if_{X_i}\mathrm{d}x_i\,\,\forall i\in\{1,2, \cdots, d\}$ and $f_{X_i}(x_i)$ is the marginal probability density function of $X_i$.
\end{definition}


\begin{definition}{Variance}{Var}
If the random variable has finite expectation $\mathbb{E}(X)$, then the expectation of $(X-\mathbb{E}(X))^2$  is called the variance of $X$ (denoted as $Var(X)$), i.e. $\sum_{i=1}^{\infty}(x_i-\mathbb{E}(X))^2P(X=x_i)$ for discrete random variable and  $\int_{-\infty}^{\infty}(x-\mathbb{E}(X))^2f_X(x)\mathrm{d}x$ for  continuous random variable.
\end{definition}

It is obvious that
$$Var(X)=\int_{\mathbb{R}}(x-\mathbb{E}(X))^2f_X(x)\mathrm{d}x=\int_{\mathbb{X}}(x-\mathbb{E}(X))^2f_X(x)\mathrm{d}x$$
where the set $\mathbb{X}$ is  the support of $X$.


In general, let $X$ be a random variable with pdf $f_X(x)$, the expectation of the random variable $g(X)$ is equal to the integration $\int_{\mathbb{R}}g(x)f_X(x)\mathrm{d}x$ if it exists.

In analysis, the expectation of a random variable is the integration of some functions.

%**********************************************************
\begin{definition}{Conditional Expectation}{Conditional Expectation}
The conditional expectation of $Y$ given $X$ can be defined as
$$\Psi(X) =\mathbb{E}(Y|X)=\int_{\mathbb{R}} y f_{Y|X}(y|x) \mathrm{d}y.$$
\end{definition}

It is a \href{https://www.encyclopediaofmath.org/index.php/Parameter-dependent_integral}{parameter-dependent integral}, where the parameter $x$ can be considered as fixed constant. See \href{http://ocw.uc3m.es/matematicas/calculus-ii/C2/unit3.pdf}{the content in Calculus II}.

\begin{definition}{Tower Rule}{Tower Rule}
\href{https://www.wikiwand.com/en/Law_of_total_expectation}{Tower Rule}
Let $X$ be random variable on a sample space $\Omega$, let the events $A_1, A_2, \cdots$ are disjoint in the sample space $\Omega$ and $\bigcup_{i=1}^{\infty}A_i = \Omega$, $B$ is any subset of $\Omega$, $\mathbb{E}(X)=\sum_{i=1}^{\infty}\mathbb{E}(X|A_i)P(A_i)$. In general, the conditional expectation $\Psi(X)=\mathbb{E}(Y|X)$ satisfies $\mathbb{E}(\Psi(X))=\mathbb{E}(Y).$
\end{definition}

The probability of an event can be expressed as the expectation, i.e.

\begin{table}[h]
	\centering
	\caption{Probability as Integration}
	\begin{tabular}{|c|}
		\hline
		$P(x\in A)= \int_{\mathbb{R}}\mathbb{I}_{A}f_X(x)\mathrm{d}x$ \\
		\hline
	\end{tabular}
\end{table}
where
\begin{equation}
\mathbb{I}_{A}=
\begin{cases}
	1, & \text{if $x\in A$} \\
	0, & \text{otherwise}
\end{cases}.
\end{equation}
It is called as the indictor function of $A$.

%******************************************************
\begin{definition}{Entropy}{Ent}
The Shannon entropy is defined as
 $$ H=-\sum_{i=1}^{\infty}P(X=x_i){\ln}P(X=x_i) $$
for the discrete random variable $X$ with the range $\{x_1,\cdots, x_n\}$ with the convention $P(X=x)=0$ that $P(X=x){\ln}{\frac{1}{P(X=x)}}=0$ and
 $$ H=\int_{\mathbb{X}}{\ln}(\frac{1}{f_X(x)})f_X(x)\mathrm{d}x $$
for the continuous random variable $X$ with the support $\mathbb{X}$.
\end{definition}

The Shannon entropy is often related to $\color{lime}{finite}$ state discrete random variable, i.e. $\color{teal}{n< \infty}$ and the value of random variable is not involved in its entropy.

\begin{table}[h]
	\centering
	\caption{Claude Elwood Shannon, 1916-2001}
	\begin{tabular}{|c|}
		\hline
		\includegraphics[width = 1.0\textwidth]{shannon.jpg}  \\
		\hline
		There is a \href{(http://colah.github.io/posts/2015-09-Visual-Information/}{Visual information} \\
        at \url{http://colah.github.io/posts/2015-09-Visual-Information/}. \\
		\hline
	\end{tabular}
\end{table}

\subsection{Moment and Moment Generating Function}

In calculus, we know that some proper functions can be extended as a Taylor series in a interval.

Moment  is a specific quantitative measure of the shape of a function.
\begin{definition}{Moments}{M}
The $n$-th moment of a real-valued continuous function $f_X(x)$ of a real variable about a value $c$ is
$$\mu_{n} =\int_{x\in\mathbb{R}}(x-c)^n f_X(x)\mathbb{d}x.$$
\end{definition}

And the constant $c$ always take the expectation $\mathbb{E}(X)$ or $0$.

\begin{definition}{Moment Generating Function}{MGF}
Moment generating function of a random variable $X$ is the expectation of the random variable of $e^{tX}$，i.e. $$M_X(t)=\mathbb{E}(e^{tX})$$.
\end{definition}

\begin{table}[h]
	\centering
	\caption{Moment generating function}
	\begin{tabular}{|c|}
		\hline
		$M_X(t)=\int_{\mathbb{R}}e^{tx}f_X(x)\mathbb{d}x$ for continuous random variable  \\
		\hline
		$M_{X}(t)=\sum_{i=1}^{\infty}e^{tx_{i}}P(x=x_i)$ for discrete random variable  \\
		\hline
	\end{tabular}
\end{table}

It is \href{http://mathworld.wolfram.com/LaplaceTransform.html}{Laplace transformation} applied to probability density function. And it is also related with cumulants \url{http://scholarpedia.org/article/Cumulants}.

The moments are not unique, i.e. the different distribution may have the same moments such as
$${f}_1(x)=\frac{1}{\sqrt{2\pi}x}\exp(-{\log(x)}^2/2)$$
$${f}_2(x)={f}_1(x)(1+\sin(2\pi\log{x})).$$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Sampling Methods}

%There are some related links in web:
%***
%
%* [**Random number generation**](https://www.wikiwand.com/en/Random_number_generation).
%* [**List of Random number generating**](https://www.wikiwand.com/en/List_of_random_number_generators).
%* [**Pseudorandom number generator**](https://www.wikiwand.com/en/Pseudorandom_number_generator).
%* [**Non-uniform pseudo-random variate generation**](https://www.wikiwand.com/en/Pseudo-random_number_sampling#/Finite_discrete_distributions).
%* [**Sampling methods**](https://ermongroup.github.io/cs228-notes/inference/sampling/).
%* [**MCMC**](https://www.wikiwand.com/en/Monte_Carlo_method).
%* [**Visualizing MCMC**](https://twiecki.github.io/blog/2014/01/02/visualizing-mcmc/).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Sampling from Discrete Distribution}

\begin{definition}{Inverse transform technique}{inverse tech}
Let $F$ be a probability cumulative function of a random variable taking non-negative integer values, and let $U$ be uniformly distributed on interval $[0, 1]$.
The random variable given by $X=k$ if and only if $F(k-1)<U<F(k)$ has the distribution function $F$.
\end{definition}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Sampling from Categorical Distribution}

The \emph{categorical distribution} (also called  a \emph{generalized Bernoulli distribution}, \emph{multinoulli distribution}) is a discrete probability distribution that describes the possible results of a random variable that can take on one of \emph{K} possible categories, with the probability of each category separately specified.

%***
The category can be represented as \href{https://www.wikiwand.com/en/One-hot#/overview}{one-hot vector}, i.e. the vector $(1,0,\cdots, 0)$ is referred as the first category and the others are as similar as it. This representation is of some advantages:
\begin{enumerate}
  \item Each category is identified by its position and equal to each other in norm;
  \item The one-hot vectors can not be compared as the real number  in value;
  \item The probability mass function of the $K$ categories distribution can be written in a compact form - $P(\mathrm{X})=[p_1,p_2, \cdots, p_K]\cdot \mathrm{X}^T$,
  where
  \begin{itemize}
   \item The probability $p_i \geq 0, \forall i\in [1,\cdots, K]$ and $\sum_{i=1}^{K}p_i=1$;
   \item The random variable $\mathrm{X}$ is one-hot vector and $\mathrm{X}^T$ is the transpose of $X$.
  \end{itemize}
\end{enumerate}
%***

Sampling it by  \href{https://www.wikiwand.com/en/Inverse_transform_sampling}{inverse transform sampling}:
\begin{enumerate}
  \item Pick a \href{https://www.wikiwand.com/en/Uniform_distribution_(continuous)}{uniformly distributed} number between 0 and 1.
  \item Locate the greatest number in the CDF whose value is less than or equal to the number just chosen., i.e. $F^{-1}(u)=inf\{x: F(x) \leq u\}$.
  \item Return the category corresponding to this CDF value.
\end{enumerate}


See more on \href{https://www.wikiwand.com/en/Categorical_distribution}{Categorical distribution}.

$\color{olive}{\text{The categorical variable cannot be ordered, how to compute the CDF?}}$

\subsection{Sampling from Continuous Distribution}

\subsubsection{Direct Methods}

\href{https://www.wikiwand.com/en/Inverse_transform_sampling)}{Inverse transform technique}
\begin{theorem}{Inverse transform technique}
Let $F$ be a probability cumulative function of a continuous random variable,
and let $U$ be uniformly distributed on interval $[0, 1]$.
The random variable $X=F^{-1}(U)$ has the distribution function $F$.
\end{theorem}

\emph{Khintchine’s (1938) theorem}: Suppose the random variable $X$ with density given by
$$ f_{X}(x)=\int_{x}^{\infty}z^{-1}f_Z(z)\mathrm{d}z $$
or This mixture can be represented equivalently by
\begin{align}
   Z &\sim f_Z(z), z > 0  &\text{and} \tag 1\\
   X|(Z=z) & \sim Uinf(0,z). \tag 2
\end{align}

Hence $\frac{X}{Z}|(Z=z)=\frac{X}{z}|(Z=z)\sim\,Unif(0,1)$ not depending on $z$, so that
\begin{align}
    \frac{X}{Z}|(Z=z) =\frac{X}{z}|(Z=z) &\sim\,Unif(0,1) \\
    \frac{X}{Z}\stackrel{d}=U &\sim Unif(0,1) \\
    X & =ZU \tag 3
\end{align}

and $U$ and $Z$ are mutually independent.

If $\nabla f_X(x)$ exists and $f_X(\infty)=0$,  we can obtain the following equation by Newton-Leibiniz's theorem
$$f_X(x)=-\int_{x}\nabla f_X(z)\mathrm{d}z $$
and comparing (4) and (0), it is obvious: $\nabla f_X(z)=-z^{-1}f_Z(z)$.

See more information on \href{https://www.jstor.org/stable/3087351?seq=1#page_scan_tab_contents}{On Khintchine's Theorem and Its Place in Random Variate Generation} and
\href{http://rspa.royalsocietypublishing.org/content/466/2119/2079}{Reciprocal symmetry, unimodality and Khintchine’s theorem}.
%%%%%#### Rejection sampling

\subsubsection{Rejection sampling}
The \href{https://www.wikiwand.com/en/Rejection_sampling}{rejection sampling} is to draw samples from a distribution with the help of a proposal distribution.
%***

%The algorithm (used by [John von Neumann](https://www.wikiwand.com/en/John_von_Neumann) and dating back to Buffon and [his needle](https://www.wikiwand.com/en/Buffon%27s_needle "Buffon's needle")) to obtain a sample from distribution $X$ with density $f$ using samples from distribution $Y$ with density $g$ is as follows:

1.   Obtain a sample $y$ from distribution $Y$ and a sample $u$ from $\mathrm{Unif}(0,1)$(the uniform distribution over the unit interval).
2.   Check whether or not $u<f(y)/Mg(y)$.
	\begin{itemize}
      \item   If this holds, accept $y$ as a sample drawn from $f$;
	  \item   if not, reject the value of $y$ and return to the sampling step.
	\end{itemize}
%***

The algorithm will take an average of $M$ iterations to obtain a sample.
\href{https://www.wikiwand.com/en/Ziggurat_algorithm}{Ziggural algorithm} is an application of rejection sampling to draw samples from Gaussian distribution.

%|[John von Neumann, 1903-1957](https://www.wikiwand.com/en/John_von_Neumann)|
%|:-------------------------------------------------------------------------:|
%|![](https://www.azquotes.com/picture-quotes/quote-anyone-who-attempts-to-generate-random-numbers-by-deterministic-means-is-of-course-living-john-von-neumann-21-31-61.jpg)|
		
\includegraphics[width = 0.9\textwidth]{john-von-neumann.jpg}  

\subsubsection{Sampling-importance-resampling}

It is also called SIR method in short. As name shown, it consists of two steps: sampling and importance sampling.
The SIR method generates an approximate i.i.d. sample of size $m$ from the target density
 $f(x),x\in \mathbb{X}\subset\mathbb{R}^n.$.
%***

The SIR without replacement is as following:
\begin{enumerate}
 	\item Draw $X^{(1)},X^{(2)},\cdots,X^{(J)}$ independently from the proposal density $g(\cdot)$.
 	\item Select a subset $\{X^{(k_i)}\}_{i=1}^{m}$ from $\{X^{(i)}\}_{i=1}^{J}$ via resampling without replacement from the discrete distribution $\{\omega_i\}_{i=1}^{J}$
    where $w_i=\frac{f(X^{(i)})}{g(X^{(i)})}$ and $\omega_i=\frac{w_i}{\sum_{i=1}^{n}w_i}$.
\end{enumerate}

\subsubsection{The Conditional Sampling Method}

The conditional sampling method due to the prominent Rosenblatt transformation is particularly available
when the joint distribution of a $d$ vector is very difficult to generate but one marginal distribution and $d-1$ univariate conditional distributions are easy to simulate.

It is based on the chain rule of probability:
$$ f(x)=f(x_d)\prod_{i=1}^{d-1}f_{k}( x_k |x_{k+1}, x_{k+2},\cdots, x_d ) $$
where $x=(x_1,x_2,\cdots, x_d)^{T}$ and $f(x)$ is the probability density function.

To generate $X$ from $f(x)$, we only need to generate $x_d$ from the marginal
density $f_d(x_d)$, then to generate $x_k$ sequentially from the conditional density
$f_k(x_k|x_{k+1}, x_{k+2}, \cdots, x_d)$.

%##### The Vertical Density Representation Method

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Numerical Optimization}

IN \href{https://homes.cs.washington.edu/~pedrod/papers/cacm12.pdf}{A Few Useful Things to Know about Machine Learning}, Pedro Domingos put up a relation:
$$\color{red}{LEARNING = REPRESENTATION + EVALUATION + OPTIMIZATION}.$$

\begin{itemize}
\item Representation as the core of the note is the general (mathematical) **model** that computer can handle.
\item Evaluation is  \textbf{criteria}. An evaluation function (also called objective function, cost function or scoring function) is needed to distinguish good classifiers from bad ones.
\item Optimization is to aimed to find the parameters that optimizes the evaluation function, i.e.
    $$ \arg\min_{\theta} f(\theta)=\{\theta^*|f(\theta^*)=\min f(\theta)\}\,\text{or}\,\arg\max_{\theta}f(\theta)=\{\theta^*|f(\theta^*)=\max f(\theta)\}. $$
\end{itemize}
%***********************************************

The objective function to be minimized is also called cost function.

Evaluation is always attached with optimization; the evaluation which cannot be optimized is not a good evaluation in machine learning.


\section{Gradient Descent and More}

Each iteration of a line search method computes a search direction $p^{k}$ and then decides how
far to move along that direction. The iteration is given by
$$ x^{k+1}=x^{k}+\alpha_{k}p^{k}\tag{Line search} $$
where the positive scalar $\alpha^{k}$ is called the step length. The success of a line search method
depends on effective choices of both the direction $p^{k}$ and the step length $\alpha_{k}$.

$\color{lime}{Note}$: we use the notation $x^{k}$ and $\alpha_k$ to represent the $k$th iteration of the vector variables $x$ and $k$th step length, respectively.
Most line search algorithms require pk to be a descent direction — one for which
$\left< {p^k},\nabla f_k \right> < 0$ — because this property guarantees that the function $f$ can be reduced along this direction, where $\nabla f_k$ is the gradient of objective function $f$ at the $k$th iteration point $x_k$ i.e. $\nabla f_k=\nabla f(x^{k})$.
%***

Gradient descent and its variants are to find the local solution of  the unconstrained optimization problem:
$$
\min f(x)
$$
where $x\in \mathbb{R}^{n}$.

Its iterative procedure is:
$$x^{k+1}=x^{k}-\alpha_{k}\nabla_{x}f(x^k)$$
where $x^{k}$ is the $k$th iterative result, $\alpha_{k}\in\{\alpha|f(x^{k+1})< f(x^{k})\}$ and particularly $\alpha_{k}=\arg\min_{\alpha}\{f(x^{k}-\alpha\nabla_{x}f(x^{k}))\}$.
%****

Some variants of gradient descent methods are not line search method.
For example, the \textbf{heavy ball method}:
$$
x^{k+1}=x^{k}-\alpha_{k}\nabla_{x}f(x^k)+\rho_{k}(x^k-x^{k-1})
$$
where the momentum coefficient $\rho_k\in[0,1]$ generally and the step length $\alpha_k$ cannot be determined by line search.

\textbf{Nesterov accelerated gradient method} is defined by

\begin{align}
x^{k}=y^{k}-\alpha^{k+1}\nabla_{x}f(y^k) \qquad &\text{Descent} \\
y^{k+1}=x^{k}+\rho^{k}(x^{k}-x^{k-1})   \qquad  &\text{Momentum}
\end{align}
where the momentum coefficient $\rho_k\in[0,1]$ generally.
\begin{table}[h]
	\centering
	\caption{Inventor of Nesterov accelerated Gradient}
	\begin{tabular}{|c|}
		\hline
		\includegraphics[width = 0.9\textwidth]{Nesterov_yurii.jpg}  \\
		\hline
	\end{tabular}
\end{table}


%* https://www.wikiwand.com/en/Gradient_descent
%* http://wiki.fast.ai/index.php/Gradient_Descent
%* https://blogs.princeton.edu/imabandit/2013/04/01/acceleratedgradientdescent/
%* https://blogs.princeton.edu/imabandit/2015/06/30/revisiting-nesterovs-acceleration/
%* http://awibisono.github.io/2016/06/20/accelerated-gradient-descent.html
%* https://jlmelville.github.io/mize/nesterov.html
%* https://smartech.gatech.edu/handle/1853/60525
%* https://zhuanlan.zhihu.com/p/41263068
%* https://zhuanlan.zhihu.com/p/35692553
%* https://zhuanlan.zhihu.com/p/35323828

\section{Mirror Gradient Method}

It is often called \textbf{mirror descent}.
It can be regarded as non-Euclidean generalization of \textbf{projected gradient descent} to solve some constrained optimization problems.

\subsection{Projected Gradient Descent}

\textbf{Projected gradient descent} is aimed to solve convex optimization problem with explicit constraints, i.e.
$$
\arg\min_{x\in\mathbb{S}}f(x)
$$
where $\mathbb{S}\subset\mathbb{R}^n$.
It has two steps:

\begin{align}
z^{k+1} = x^{k}-\alpha_k\nabla_x f(x^{k}) &\qquad \text{Gradient descent}\\
x^{k+1} = Proj_{\mathbb{S}}(z^{k+1})=\arg\min_{x\in \mathbb{S}}\|x-z^{k+1}\|^{2} &\qquad\text{Projection}
\end{align}


\subsection{Mirror Descent}

\textbf{Mirror descent} can be regarded as the non-Euclidean generalization via replacing the $\ell_2$ norm or Euclidean distance in projected gradient descent by \href{https://www.mdpi.com/1099-4300/16/12/6338/htm}{Bregman divergence}.

Bregman divergence is induced by convex smooth function $f$:
$$
B(x,y)=f(x)-f(y)-\left<\nabla f(y),x-y\right>
$$
where $\left<\cdot,\cdot\right>$ is inner product.
Especially, when $f$ is quadratic function, the Bregman divergence induced by $f$ is
$$
B(x,y)=x^2-y^2-\left<2y,x-y\right>=x^2+y^2-2xy=(x-y)^2
$$
i.e. the Euclidean distance.
A wonderful introduction to Bregman divergence is \textbf{Meet the Bregman Divergences} by \href{http://mark.reid.name}{Mark Reid} at \url{http://mark.reid.name/blog/meet-the-bregman-divergences.html}.
%***
It is given by:
\begin{align}
z^{k+1} = x^{k}-\alpha_k\nabla_x f(x^{k}) &\qquad \text{Gradient descent}\\
x^{k+1} = \arg\min_{x\in\mathbb{S}}B(x,z^{k+1}) &\qquad\text{Bregman projection}
\end{align}

One special method is called \textbf{entropic mirror descent} when $f=e^x$ and $\mathbb{S}$ is simplex.

See more on the following link list.

%* http://users.cecs.anu.edu.au/~xzhang/teaching/bregman.pdf
%* https://zhuanlan.zhihu.com/p/34299990
%* https://blogs.princeton.edu/imabandit/2013/04/16/orf523-mirror-descent-part-iii/
%* https://blogs.princeton.edu/imabandit/2013/04/18/orf523-mirror-descent-part-iiii/
%* https://www.stat.berkeley.edu/~bartlett/courses/2014fall-cs294stat260/lectures/mirror-descent-notes.pdf

\section{Variable Metric Methods}

\subsection{Newton's Method}

NEWTON’S METHOD and QUASI-NEWTON METHODS are classified to variable metric methods.

It is also to find the solution of unconstrained optimization problems, i.e.
$$\min f(x)$$
where $x\in \mathbb{R}^{n}$.
%***

\textbf{Newton's method} is given by
$$
x^{k+1}=x^{k}-\alpha^{k+1}H^{-1}(x^{k})\nabla_{x}\,{f(x^{k})}
$$
where $H^{-1}(x^{k})$ is inverse of the Hessian matrix of the function $f(x)$ at the point $x^{k}$.
It is called \textbf{Newton–Raphson algorithm} in statistics.
Especially when the log-likelihood function $\ell(\theta)$ is well-behaved,
a natural candidate for finding the MLE is the Newton–Raphson algorithm with quadratic convergence rate.


\subsection{The Fisher Scoring Algorithm}

In maximum likelihood estimation, the objective function is the log-likelihood function, i.e.
$$
\ell(\theta)=\sum_{i=1}^{n}\log{P(x_i|\theta)}
$$
where $P(x_i|\theta)$ is the probability of realization $X_i=x_i$ with the unknown parameter $\theta$.
However, when the sample random variable $\{X_i\}_{i=1}^{n}$ are not observed or realized, it is best to
replace negative Hessian matrix (i.e. -$\frac{\partial^2\ell(\theta)}{\partial\theta\partial\theta^{T}}$) of the likelihood function with the
\textbf{observed information matrix}:
$$
J(\theta)=\mathbb{E}(\color{red}{\text{-}}\frac{\partial^2\ell(\theta)}{\partial\theta\partial\theta^{T}})=\color{red}{-}\int\frac{\partial^2\ell(\theta)}{\partial\theta\partial\theta^{T}}f(x_1, \cdots, x_n|\theta)\mathrm{d}x_1\cdots\mathrm{d}x_n
$$
where $f(x_1, \cdots, x_n|\theta)$ is the joint probability density function of  $X_1, \cdots, X_n$ with unknown parameter $\theta$.

And the \textbf{Fisher scoring algorithm} is given by
$$\theta^{k+1}=\theta^{k}+\alpha_{k}J^{-1}(\theta^{k})\nabla_{\theta} \ell(\theta^{k})$$
where $J^{-1}(\theta^{k})$ is the inverse of observed information matrix at the point $\theta^{k}$.

See \url{http://www.stats.ox.ac.uk/~steffen/teaching/bs2HT9/scoring.pdf} or \url{https://wiseodd.github.io/techblog/2018/03/11/fisher-information/} for more information.

Fisher scoring algorithm is regarded  as an example of \textbf{Natural Gradient Descent} in
information geometry  such as \url{https://wiseodd.github.io/techblog/2018/03/14/natural-gradient/}.

\subsection{Quasi-Newton Methods}

Quasi-Newton methods, like steepest descent, require only the gradient of the objective
function to be supplied at each iterate.
By measuring the changes in gradients, they construct a model of the objective function
that is good enough to produce superlinear convergence.
The improvement over steepest descent is dramatic, especially on difficult
problems. Moreover, since second derivatives are not required, quasi-Newton methods are
sometimes more efficient than Newton's method.

In optimization, quasi-Newton methods (a special case of 、textbf{variable-metric methods}) are algorithms for finding local maxima and minima of functions. 
Quasi-Newton methods are based on Newton's method to find the stationary point of a function, where the gradient is $0$.
In quasi-Newton methods the Hessian matrix does not need to be computed. The Hessian is updated by analyzing successive gradient vectors instead. Quasi-Newton methods are a generalization of the secant method to find the root of the first derivative for multidimensional problems. In multiple dimensions the secant equation is under-determined, and quasi-Newton methods differ in how they constrain the solution, typically by adding a simple low-rank update to the current estimate of the Hessian.
One of the chief advantages of quasi-Newton methods over Newton's method is that the Hessian matrix (or, in the case of quasi-Newton methods, its approximation) $B$ does not need to be inverted. The Hessian approximation $B$ is chosen to satisfy
$$ \nabla f(x^{k+1})=\nabla f(x^{k})+B(x^{k+1}-x^{k}), $$
which is called the \textbf{secant equation} (the Taylor series of the gradient itself).
In more than one dimension $B$ is underdetermined. In one dimension, solving for B and applying the Newton's step with the updated value is equivalent to the \href{https://www.wikiwand.com/en/Secant_method}{secant method}.
The various quasi-Newton methods differ in their choice of the solution to the \emph{secant equation} (in one dimension, all the variants are equivalent).

%![BFGS](http://aria42.com/images/bfgs.png)
%***
%
%* [Wikipedia page](https://www.wikiwand.com/en/Newton%27s_method_in_optimization)
%* [Newton-Raphson Visualization (1D)](http://bl.ocks.org/dannyko/ffe9653768cb80dfc0da)
%* [Newton-Raphson Visualization (2D)](http://bl.ocks.org/dannyko/0956c361a6ce22362867)
%* [Newton's method](https://www.wikiwand.com/en/Newton%27s_method)
%* [Quasi-Newton method](https://www.wikiwand.com/en/Quasi-Newton_method)
%* [Using Gradient Descent for Optimization and Learning](http://www.gatsby.ucl.ac.uk/teaching/courses/ml2-2008/graddescent.pdf)
%* http://fa.bianp.net/teaching/2018/eecs227at/quasi_newton.html

\subsection{Natural Gradient Descent}

\textbf{Natural gradient descent} is to solve the optimization problem $\min_{\theta} L(\theta)$ by
$$ \theta^{(t+1)}=\theta^{(t+1)}-\alpha_{(t)}F^{-1}(\theta^{(t)})\nabla_{\theta}L(\theta^{(t)}) $$
where $F^{-1}(\theta^{(t)})$ is the inverse of Remiann metric at the point $\theta^{(t)}$.
And Fisher scoring algorithm is a typical application of Natural Gradient Descent to statistics.
\textbf{Natural gradient descent} for manifolds corresponding to
exponential families can be implemented as a first-order method through \emph{mirror descent} at \url{https://www.stat.wisc.edu/~raskutti/publication/MirrorDescent.pdf}.

\begin{table}[h]
	\centering
	\caption{Originator of Information Geometry: Shun-ichi Amari(甘利 俊一)}
	\begin{tabular}{|c|}
		\hline
		\includegraphics[width = 0.80\textwidth]{shun-ichi-amari.jpg}  \\
		\hline
	\end{tabular}
\end{table}
%| Originator of Information Geometry |
%|:----:|
%|![Shun-ichi Amari](https://groups.oist.jp/sites/default/files/imce/u34/images/people/shun-ichi-amari.jpg)|

%* http://www.yann-ollivier.org/rech/publs/natkal.pdf
%* http://www.dianacai.com/blog/2018/02/16/natural-gradients-mirror-descent/
%* https://www.zhihu.com/question/266846405
%* http://bicmr.pku.edu.cn/~dongbin/Conferences/Mini-Course-IG/index.html
%* http://ipvs.informatik.uni-stuttgart.de/mlr/wp-content/uploads/2015/01/mathematics_for_intelligent_systems_lecture12_notes_I.pdf
%* http://www.luigimalago.it/tutorials/algebraicstatistics2015tutorial.pdf
%* http://www.yann-ollivier.org/rech/publs/tango.pdf
%* http://www.brain.riken.jp/asset/img/researchers/cv/s_amari.pdf

\section{ Expectation Maximization Algorithm}

\textbf{Expectation-Maximization algorithm}, popularly known as the  EM algorithm has become a standard piece in the statistician's repertoire.
It is used in incomplete-data problems or latent-variable problems such as Gaussian mixture model in maximum likelihood  estimation.
The basic principle behind the EM is that instead of performing a complicated optimization,
one augments the observed data with latent data to perform a series of simple optimizations.

Let $\ell(\theta|Y_{obs})\stackrel{\triangle}=\log{L(\theta|Y_{obs})}$ denote the log-likelihood function of observed datum $Y_{obs}$.
We augment the observed data $Y_{obs}$ with latent variables $Z$ so that both the
complete-data log-likelihood $\ell(\theta|Y_{obs}, Z)$ and the conditional predictive distribution $f(z|Y_{obs}, \theta)$ are available.
Each iteration of the EM algorithm consists of an expectation step (E-step) and a maximization step (M-step)
Specifically, let $\theta^{(t)}$ be the current best guess at the MLE $\hat\theta$. The E-step
is to compute the \textbf{Q} function defined by

\begin{align}
Q(\theta|\theta^{(t)})
&= \mathbb{E}(\ell(\theta|Y_{obs}, Z)|Y_{obs},\theta^{(t)}) \\
&= \int_{Z}\ell(\theta|Y_{obs}, Z)\times f(z|Y_{obs}, \theta^{(t)})\mathrm{d}z,
\end{align}
and the M-step is to maximize **Q** with respect to $\theta$ to obtain
$$\theta^{(t+1)}=\arg\max_{\theta} Q(\theta|\theta^{(t)}).$$

%* https://www.wikiwand.com/en/Expectation%E2%80%93maximization_algorithm
%* http://cs229.stanford.edu/notes/cs229-notes8.pdf
%* https://www2.stat.duke.edu/courses/Spring06/sta376/Support/EM/EM.Mixtures.Figueiredo.2004.pdf

\begin{table}[h]
	\centering
	\caption{Diagram of EM algorithm}
	\begin{tabular}{|c|}
		\hline
		\includegraphics[width = 1.0\textwidth]{EM.png}  \\
		\hline
	\end{tabular}
\end{table}
%|Diagram of EM algorithm|
%|:---------------------:|
%|![](https://i.stack.imgur.com/v5bqe.png)|
%EM.png

\subsection{Generalized EM Algorithm}

Each iteration of the \textbf{generalized EM} algorithm consists of an expectation step (E-step) and an ascent step instead of maximization step (M-step).
Specifically, let $\theta^{(t)}$ be the current best guess at the MLE $\hat\theta$. The E-step
is to compute the **Q** function defined by
\begin{align}
Q(\theta|\theta^{(t)}) &= \mathbb{E}(\ell(\theta|Y_{obs}, Z)|Y_{obs},\theta^{(t)}) \\
                       &= \int_{Z}\ell(\theta|Y_{obs}, Z)\times f(z|Y_{obs}, \theta^{(t)})\mathrm{d}z,
\end{align}

and the another step is to find  $\theta$ that satisfies $Q(\theta^{t+1}|\theta^{t})>Q(\theta^{t}|\theta^{t})$, i.e.
$$\theta^{(t+1)}\in \{\hat{\theta}|Q(\hat{\theta}|\theta^{(t)} \geq Q(\theta|\theta^{(t)}) \}.$$

It is not to maximize the conditional expectation.

See more on the book \textbf{The EM Algorithm and Extensions, 2nd Edition
by Geoffrey McLachlan , Thriyambakam Krishna} at \url{https://www.wiley.com/en-cn/The+EM+Algorithm+and+Extensions,+2nd+Edition-p-9780471201700}.

See more on \url{https://www.stat.berkeley.edu/~aldous/Colloq/lange-talk.pdf}

\section{Alternating Direction Method of Multipliers}

Alternating direction method of multipliers is called \textbf{ADMM} shortly.
It is aimed to solve the following convex optimization problem:
\begin{align}
  \min F(x,y) \{&=f(x)+g(y)\} \tag {cost function} \\
          Ax+By &=b           \tag{constraint}
\end{align}

where $f(x)$ and $g(y)$ are convex; $A$ and $B$ are matrices.

Define the augmented Lagrangian:
$$ L_{\mu}(x,y)=f(x)+g(y)+\lambda^{T}(Ax+By-b)+\frac{\mu}{2}\|Ax+By-b\|_{2}^{2}. $$
%***

It is iterative procedure at $k$th step:
\begin{enumerate}
\item $x^{k+1}=\arg\min_{x}L_{\mu}(x,y^{\color{teal}{k}},\lambda^{\color{teal}{k}});$
\item $y^{k+1}=\arg\min_{y} L_{\mu}(x^{\color{red}{k+1}}, y, \lambda^{\color{teal}{k}});$
\item $\lambda^{k+1} = \lambda^{k}+\mu (Ax^{\color{red}{k+1}} + By^{\color{red}{k+1}}-b).$
\end{enumerate}

\begin{algorithm}
	\caption{ADMM}
	\label{ADMM}
	\begin{algorithmic}
	\REQUIRE $\text{the initial points:}\quad x_0, y_0, {\lambda}_0$
	\ENSURE $k = 0$
	\WHILE{$k \neq N$}
	%\IF{$N$ is even}
	\STATE $x^{k+1}=\arg\min_{x}L_{\mu}(x,y^{\color{teal}{k}},\lambda^{\color{teal}{k}})$
	\STATE $y^{k+1}=\arg\min_{y} L_{\mu}(x^{\color{red}{k+1}}, y, \lambda^{\color{teal}{k}})$
	\STATE $\lambda^{k+1} = \lambda^{k}+\mu (Ax^{\color{red}{k+1}} + By^{\color{red}{k+1}}-b)$
	\STATE $k \Leftarrow  k + 1$
	%\ENDIF
	\ENDWHILE
	\end{algorithmic}
\end{algorithm}


%$\color{teal}{\text{Thanks to Professor He Bingsheng who taught me this.}}$
\begin{note}
	Thanks to \href{http://maths.nju.edu.cn/~hebma/}{Professor He Bingsheng} who taught me this.
\end{note}
%* https://www.ece.rice.edu/~tag7/Tom_Goldstein/Split_Bregman.html
%* http://maths.nju.edu.cn/~hebma/
%* http://stanford.edu/~boyd/admm.html
%* http://shijun.wang/2016/01/19/admm-for-distributed-statistical-learning/
%* https://www.wikiwand.com/en/Augmented_Lagrangian_method
%* https://blog.csdn.net/shanglianlm/article/details/45919679
%* http://www.optimization-online.org/DB_FILE/2015/05/4925.pdf


\section{ Stochastic Gradient Descent}

\textbf{Stochastic gradient descent} takes advantages of stochastic or estimated gradient to replace the true gradient in gradient descent.
It is  \emph{stochastic gradient} but may not be \emph{descent}.
The name \textbf{stochastic gradient methods}  may be more appropriate to call the methods with stochastic gradient.
It can date back to \textbf{stochastic approximation}.

It is aimed to solve the problem with finite sum optimization problem, i.e.
$$\arg\min_{\theta}\frac{1}{n}\sum_{i=1}^{n}f(\theta|x_i)$$
where $n<\infty$ and $\{f(\theta|x_i)\}_{i=1}^{n}$ are in the same function family and $\{x_i\}_{i=1}^{n}\subset \mathbb{R}^{d}$ are constants
while $\theta\in\mathbb{R}^{p}$ is the variable vector.

The difficulty is that $p$, that the dimension of $\theta$, is tremendous. In other words, the model is \textbf{overparameterized}.
And the number $n$ is far larger than $p$ generally, i.e. $n \gg  p\gg d$.
What is worse, the functions  $\{f(\theta|x_i)\}_{i=1}^{n}$ are not convex in most case.
%***

The stochastic gradient method is defined as
$$ \theta^{k+1}=\theta^{k}-\alpha_{k}\frac{1}{m}\sum_{j=1}^{m}\nabla f(\theta^{k}| x_{j}^{\prime}) $$
where $x_{j}^{\prime}$ is stochastically draw from $\{x_i\}_{i=1}^{n}$ and $m\ll n$.

It is the fact $m\ll n$ that makes it possible to compute the gradient of finite sum objective function and its side effect is that the objective function is not always descent.
There is fluctuations in the total objective function as gradient steps with respect to mini-batches are taken.

%******************************************************************
\begin{table}[h]
	\centering
	\caption{The fluctuations in the objective function as gradient with respect to mini-batches are taken}
	\begin{tabular}{|c|}
		\hline
		\includegraphics[width = 0.90\textwidth]{SGD.png}  \\
		\hline
	\end{tabular}
\end{table}

%|The fluctuations in the objective function as gradient steps with respect to mini-batches are taken|
%|:------------------------------------:|
%|![](https://upload.wikimedia.org/wikipedia/commons/f/f3/Stogra.png)|


An heuristic proposal for avoiding the choice and for modifying the learning rate while the learning task runs is the bold driver (BD) method.
The learning rate increases exponentially if successive steps reduce the objective function $f$, 
and decreases rapidly if an ``accident” is encountered (if objective function $f$ increases), until a suitable value is found.
After starting with a small learning rate, its modifications are described by the following equation:
$$
\alpha_{k+1}=
   \begin{cases}
       \rho \alpha_{k}, & {f(\theta^{k+1})< f(\theta^{k})}; \\
       \eta^n \alpha_{k}, & {f(\theta^{k+1})> f(\theta^{k})} \text{using ${\alpha}_k$},
   \end{cases}
$$
where $\rho$ is close to 1 such as $\rho=1.1$  in order to avoid frequent “accidents” because the
objective function computation is wasted in these cases, $\eta$ is chosen to provide a rapid reduction
($\eta = 0.5$), and $n$ is the minimum integer such that the reduced rate $\eta^n$ succeeds in diminishing the objective function.
%***

The fact that the sample size is far larger than the dimension of parameter, $n\gg p$,
that makes it expensive to compute total objective function $f(\theta)=\sum_{i=1}^{n}f(\theta|{x}_i)$.
Thus it is not clever to determine the learning rate $\alpha_k$ by line search.
And most stochastic gradient methods are to find  proper step length $\alpha_{k}$ to make it converge at least in convex optimization.
The variants of gradient descent such as momentum methods or mirror gradient methods have their stochastic counterparts.
\begin{itemize}
\item It is simplest to set the step length a constant, such as ${\alpha}_k=3\times 10^{-3}\, \forall k$.
\item There are decay schemes, i.e. the step length ${\alpha}_k$ diminishes such as ${\alpha}_k=\frac{\alpha}{k}$, where $\alpha$ is constant.
\item And another strategy is to tune the step length adaptively such as \textbf{AdaGrad, ADAM}.
\end{itemize}

$\color{lime}{PS}$: the step length  $\alpha_k$ is called **learning rate** in machine learning and
stochastic gradient descent is also named as \emph{incremental gradient descent method} such as the survey at \url{http://www.mit.edu/~dimitrib/Incremental_Survey_LIDS.pdf} in some case.

See the following links for more information on \emph{stochastic gradient descent}.

%\animategraphics{24}{sgd-}{0}{300}

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=0.9\textwidth]{sgdvsgd.png}
	\caption{The Differences of Gradient Descent and Stochastic Gradient Descent\label{fig:sgd vs gd}}
\end{figure}

%* https://www.wikiwand.com/en/Stochastic_gradient_descent
%* https://www.bonaccorso.eu/2017/10/03/a-brief-and-comprehensive-guide-to-stochastic-gradient-descent-algorithms/
%* https://leon.bottou.org/projects/sgd
%* https://leon.bottou.org/research/stochastic
%* http://ruder.io/optimizing-gradient-descent/
%* http://ruder.io/deep-learning-optimization-2017/
%* https://zhuanlan.zhihu.com/p/22252270
%* https://henripal.github.io/blog/stochasticdynamics
%* https://henripal.github.io/blog/langevin
%* http://fa.bianp.net/teaching/2018/eecs227at/stochastic_gradient.html
%
%|The Differences of Gradient Descent and Stochastic Gradient Descent|
%|:-----------------------------------------------------------------:|
%|![](https://wikidocs.net/images/page/3413/sgd.png)|
%
%***
%
%|The Differences of Stochastic Gradient Descent and its Variants|
%|:-------------------------------------------------------------:|
%|![](http://beamandrew.github.io//images/deep_learning_101/sgd.gif)|
%
%- [ ] http://www.optimization-online.org/
%- [ ] More [Optimization Online Links](http://www.optimization-online.org/links.html)
%- [ ] **TUTORIALS AND BOOKS** at <http://plato.asu.edu/sub/tutorials.html>.
%- [ ] http://sunju.org/research/nonconvex/
%- [ ] https://arxiv.org/pdf/1712.07897.pdf
%- [ ] https://arxiv.org/pdf/1707.02444.pdf
%- [ ] http://www.vision.jhu.edu/assets/HaeffeleCVPR17.pdf
%- [ ] https://core.ac.uk/display/73408878
%- [ ] https://www.springer.com/us/book/9783319314822
%- [ ] https://core.ac.uk/display/83849901
%- [ ] https://zhuanlan.zhihu.com/p/51514687
%- [ ] http://math.cmu.edu/~hschaeff/research.html

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Artificial Neural Network}

Artificial neural networks are most easily visualized in terms of a \textbf{directed graph}. 
In the case of sigmoidal units, node $s$ represents sigmoidal unit  and directed edge $e=(u,v)$ indicates that one of sigmoidal unit $v$'s inputs is the output of sigmoidal unit $u$.
\href{https://brilliant.org/wiki/artificial-neural-network/}{A brief introduction to artificial neural network} tells the basics of artificial neural network in a visual intuitive way.


\clearpage
\begin{figure}[!htbp]
	\centering
	\includegraphics[width=\textwidth]{neuralnetworks.png}
	\caption{Neural Networks\label{fig:nn}}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Perceptron}

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=\textwidth]{perceptron.png}
	\caption{Perceptron\label{fig:perceptron}}
\end{figure}
 
It is can be decomposed into $2$ steps:
\begin{enumerate}
	\item Aggregate all the information: $z=\sum_{i=1}^{n}w_ix_i+b_0=(x_1,x_2,\cdots,x_n,1)\cdot (w_1,w_2,\cdots,w_n,b_0)^{T}.$
	\item Transform the information to activate something: $y=\sigma(z)$, where $\sigma$ is nonlinear such as step function.
\end{enumerate}

\href{http://www.cs.columbia.edu/~mcollins/courses/6998-2012/notes/perc.converge.pdf}{Convergence Proof for the Perceptron Algorithm} shows the proof that perceptron is linear classifier.

\subsection{Activation function}

The nonlinear function $\sigma$ is conventionally called activation function.
There are some activation functions in history.
The activation function $\sigma$ can extend to more nonlinear function.
\begin{center}
\begin{itemize}
\item Sign function
   $$f(x)=
   \begin{cases}
	     1,&\text{if $x > 0$}\\
		-1,&\text{if $x < 0$}
	\end{cases}
   $$
\item Step function
	$$f(x)=
	\begin{cases}
		1,&\text{if $x\geq0$}\\
		0,&\text{otherwise}
	\end{cases}
    $$
\item Sigmoid function
    $$\sigma(x)=\frac{1}{1+e^{-x}}.$$

\item Radical base function
    $$\rho(x)=e^{-\beta(x-x_0)^2}.$$

\item TanH function
   $$tanh(x)=2\sigma(2x)-1=\frac{2}{1+e^{-2x}}-1.$$

\item ReLU function
   $$ReLU(x)={(x)}_{+}=\max\{0,x\}=
		\begin{cases}
			x,&\text{if $x\geq 0$}\\
			0,&\text{otherwise}
		\end{cases}.
   $$
\end{itemize}
\end{center}


\section{Feed-forward Neural Networks}

\subsection{Representation of Feedforward Neural Network}

Given that the function of a single neuron is rather simple, it subdivides the input space into two regions by a hyperplane, the complexity must come from having more layers of neurons involved in a complex action (like recognizing your grandmother in all possible situations).The "squashing" functions introduce critical nonlinearities in the system, without their presence multiple layers would still create linear functions.
Organized layers are very visible in the human cerebral cortex, the part of our brain which plays a key role in memory, attention, perceptual awareness, thought, language, and consciousness.%[^13]

The \textbf{feed-forward neural network} is also called multilayer perceptron. 
\href{http://math.mit.edu/~gs/learningfromdata/SIAM03.pdf}{The best way to create complex functions from simple functions is by \emph{composition}.}

In mathematics, it  can be considered as multi-layered non-linear composite function:
$$X\to \sigma\circ (W_1X+b_1)=H_1\to \sigma\circ(W_2 H_1+b_2)=H_2 \to\cdots \sigma(WH+b)=y $$
where the notation $\circ$, $M_1,b_1,M_2,b_2,\cdots, W,b$ mean pointwise operation, the parameters in the affine mapping, respectively. 
Thus the data flow in the form of the chain:
\begin{align}
f = H_1 \circ {H_2}\circ \cdots \circ{\sigma} \qquad   &\text{Composite form}           \\
X\stackrel{\sigma}\to H_1 \stackrel{\sigma}\to H_2 \stackrel{\sigma} \to \cdots\stackrel{\sigma}\to\,y \qquad &\text{Hierarchy form}         \\
\mathbb{R}^p\to \mathbb{R}^{l_1}\to \mathbb{R}^{l_2}\to \cdots\to \mathbb{R}  \qquad&  \text{Dimension}
\end{align}
where the circle notation $\circ$ means forward composite or as the input of afterward operation.
In hierarchy form, we omit the affine map.

\begin{definition}{Feed-forward Neural Networks}{fnn}
	The feedforward neural networks are written in the \emph{recursive form}:
	\begin{align}
		\mathbf{z}_{i} &= W_{i}H_{i-1}+b_i, \\
                H_{i}  &=\sigma\circ (\mathbf{z}_{i}),
		\forall i\{1,2,\dots,D\} 
	\end{align}
where $H_{0}=X\in\mathbb{R}^p$, $b_i$ is a vector and $W_{i}$ is matrix. And the number of recursive times $D$ is called the depth of network.
\end{definition}

\begin{note}
	We can compute the input in each layer and sent the output to the next layer.
	\begin{itemize}
		\item In the first layer, we feed the input vector $X\in\mathbb{R}^{p}$ and connect it to each unit in the next layer $W_1X+b_1\in\mathbb{R}^{l_1}$ 
		      where $W_1\in\mathbb{R}^{n\times l_1}, b_1\in\mathbb{R}^{l_1}$. 
              The output of the first layer is $H_1=\sigma\circ(M_1X+b)$, or in another word the output of $j$th unit in the first (hidden) layer is $h_j=\sigma{(W_1X+b_1)}_j$ 
               where ${(W_1X+b_1)}_j$ is the $j$th element of $l_1$-dimensional vector $W_1X+b_1$.
        \item In the second layer, its input is the output of first layer,$H_1$, and apply linear map to it: $W_2H_1+b_2\in\mathbb{R}^{l_2}$, 
               where $W_2\in\mathbb{R}^{l_1\times l_2}, b_2\in\mathbb{R}^{l_2}$. The output of the second layer is $H_2=\sigma\circ(W_2H_1+b_2)$, 
                or in another word the output of $j$th unit in the second (hidden) layer is $h_j=\sigma{(W_2H_1+b_2)}_j$ where ${(W_2H_1+b_2)}_j$ is the $j$th element of $l_2$-dimensional vector $W_2H_1+b_2$.
		\item The map between the second layer and the third layer is similar to (1) and (2): the linear maps datum to different dimensional space and the nonlinear maps extract better representations.
		\item In the last layer, suppose the input data is $H\in\mathbb{R}^{l}$. The output may be vector or scalar values and $W$ may be a matrix or vector as well as $y$.
	\end{itemize}
\end{note}

\par
The ordinary feedforward neural networks take the *sigmoid* function $\sigma(x)=\frac{1}{1+e^{-x}}$ as the nonlinear activation function
while the \emph{RBF networks} take the \href{https://www.wikiwand.com/en/Radial_basis_function}{Radial basis function} as the activation function such as $\sigma(x)=e^{c{\|x\|}_2^2}$.

In theory, the universal approximation theorem show the power of feedforward neural network
if we take some proper activation functions such as sigmoid function.

%* https://www.wikiwand.com/en/Universal_approximation_theorem
%* http://mcneela.github.io/machine_learning/2017/03/21/Universal-Approximation-Theorem.html
%* http://neuralnetworksanddeeplearning.com/chap4.html

\subsection{Evaluation and Optimization in Multilayer Perceptron}

The problem is how to find the optimal parameters $W_1, b_1, W_2, b_2,\cdots, W, b$ ?
The multilayer perceptron is as one example of supervised learning, which means that we feed datum
$ D=\{(\mathbf{x_i},d_i)\}_{i=1}^{n}$ to it and evaluate it.

The general form of the evaluation is given by:
$$ J(\theta)=\frac{1}{n}\sum_{i=1}^{n}\mathbb{L}[f(\mathbf{x}_i|\theta),\mathbf{d}_i]$$
where $\mathbf{d}_i$ is the desired value of the input $\mathbf{x}_i$ and $\theta$ is the parameters of multilayer perceptron. 
The notation $f(\mathbf{x}_i|\theta)$ is the output given parameters $\theta$. 
The function $\mathbb{L}$ is \emph{loss function} to measure the discrepancy between the predicted value $f(\mathbf{x}_i|\theta)$ and the desired value $\mathbf{d}_i$.

In general, the number of parameters $\theta$ is less than the sample size $n$. And the objective function $J(\theta)$ is not convex.  

We will solve it in the next section \emph{Backpropagation, Optimization and Regularization}.

%***

\href{https://arogozhnikov.github.io/3d_nn/}{Visualizing level surfaces of a neural network with raymarching}

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=\textwidth]{mlp.jpg}
	\caption{The diagram of MLP\label{fig:MLP}}
\end{figure}


%|The diagram of MLP|
%|:----------------:|
%|![](https://www.hindawi.com/journals/jcse/2012/389690.fig.0011.jpg)|
%|[Visualizing level surfaces of a neural network with raymarching](https://arogozhnikov.github.io/3d_nn/)|
%

%* https://devblogs.nvidia.com/deep-learning-nutshell-history-training/
%* [Deep Learning 101 Part 2](http://beamandrew.github.io/deeplearning/2017/02/23/deep_learning_101_part2.html)
%* https://www.wikiwand.com/en/Multilayer_perceptron
%* https://www.wikiwand.com/en/Feedforward_neural_network
%* https://www.wikiwand.com/en/Radial_basis_function_network

\subsection{Evaluation for different tasks}

Evaluation is to judge the models with different parameters in some sense via objective function.
In maximum likelihood estimation, it is likelihood or log-likelihood;
in parameter estimation, it is bias or mean square error;
in regression, it depends on the case.

In machine learning, evaluation is aimed to  measure the discrepancy between the predicted values and trues value by **loss function**.
It is expected to be continuous smooth and differential but it is snot necessary.
The principle is the loss function make the optimization or learning tractable.
For classification, the loss function always is cross entropy;
for regression, the loss function can be any norm function such as the $\ell_2$ norm;
in probability models, the loss function always is joint probability or logarithm of joint probability.

In classification, the last layer is to predict the degree of belief of the labels via \href{http://freemind.pluskid.org/machine-learning/softmax-vs-softmax-loss-numerical-stability/}{softmax function},
i.e.
$$softmax(z)=(\frac{\exp(z_1)}{\sum_{i=1}^{n}\exp(z_i)},\frac{\exp(z_2)}{\sum_{i=1}^{n} \exp(z_i)}, \cdots, \frac{\exp(z_n)}{\sum_{i=1}^{n}\exp(z_i)})$$
where $n$ is the number of total classes. 
The labels are encoded as the one hot vector such as $\mathrm{d}=(1,0,0,\cdots,0)$. The \href{https://www.cnblogs.com/smartwhite/p/8601477.html}{cross entropy} is defined as:
$$ \mathbf{H}(d,p)=-\sum_{i=1}^{n}d_i\log(p_i)=\sum_{i=1}^{n}d_i\log(\frac{1}{p_i}),$$
where $d_i$ is the $i$th element of the one-hot vector $d$ and $p_i=\frac{\exp(z_i)}{\sum_{j=1}^{n}\exp(z_j)}$ for all $i=1,2, \dots, n.$

Suppose $\mathrm{d}=(1,0,0,\cdots,0)$, the cross entropy is $\mathbf{H}(d,p)=-\log(p_1)=\log \sum_{i=1}^{n}\exp(z_i)-z_1$. 
The cost function is $\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}(d^{i},p^{i})$ in the training data set $\{(\mathbf{x}_i,d^i)\}_{i=1}^{n}$ where $\mathbf{x}_i$ is the features of $i$th sample and $d^i$ is the desired true target label encoded in **one-hot** vector meanwhile $p^{i}$ is the predicted label of $\mathbf{x}_i$.
%See the following links for more information on cross entropy and softmax.

\begin{figure}[H]
\centering
  \subfigure[SubfigureCaption]{
    \label{Fig.sub.1}
    \includegraphics[width=0.45/textwidth]{resnet56_noshort_small.jpg}}
  \subfigure[SubfigureCaption]{
    \label{Fig.sub.2}
    \includegraphics[width=0.45/textwidth]{resnet56_small.jpg}}
  \caption{VISUALIZING THE LOSS LANDSCAPE OF NEURAL NETS|}
\label{fig: visual loss}
\end{figure}

%|VISUALIZING THE LOSS LANDSCAPE OF NEURAL NETS||
%|:-------------------------------------------:|---|
%|![VGG](https://raw.githubusercontent.com/tomgoldstein/loss-landscape/master/doc/images/resnet56_noshort_small.jpg)|![ResNet](https://raw.githubusercontent.com/tomgoldstein/loss-landscape/master/doc/images/resnet56_small.jpg)|

%* <https://blog.csdn.net/u014380165/article/details/77284921>;
%* <https://blog.csdn.net/u014380165/article/details/79632950>;
%* <https://eli.thegreenplace.net/2016/the-softmax-function-and-its-derivative/>;
%* <https://www.zhihu.com/question/65288314>.

In regression, the loss function may simply be the squared $\ell_2$ norm, i.e. $\mathbb{L}(d,p)=(d-p)^{2}$ where $d$ is the desired target and $p$ is the predicted result. 
And the cost function is \emph{mean squared error}:
$$ J(\theta)=\frac{1}{n}\sum_{i=1}^{n}[f(\mathbf{x}_i|\theta)-\mathrm{d}_i]^2.$$
In \textbf{robust statistics}, there are more loss functions such as \emph{Huber loss, hinge loss, Tukey loss}.
%***
\par
\begin{itemize}
\item \href{https://www.wikiwand.com/en/Huber_loss}{Huber loss function}
   $$Huber_{\delta}(x)=
        \begin{cases}
            \frac{|x|}{2},&\text{if $|x|\leq\delta$}\\
            \delta(|x|-\frac{1}{2}\delta),&\text{otherwise}
        \end{cases}
   $$

\item \href{https://www.wikiwand.com/en/Hinge_loss}{Hinge loss function}
  $$Hinge(x)=max\{0, 1-tx\}$$
  where $t=+1$ or $t=-1$.

\item Tukey loss function
   $$
	Tukey_{\delta}(x)=
	\begin{cases}
         (1-[1-x^2/\delta^2]^3)\frac{\delta^2}{6},&\text{if $|x|\leq\delta$}\\
         \frac{\delta^2}{6},                      &\text{otherwise}
    \end{cases}
   $$

\end{itemize}
\par

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=0.6\textwidth]{RhoFunctions.png}
	\caption{LOSS FUNCTION\label{fig:loss}}
\end{figure}

It is important to choose or design loss function or more generally objective function,
which can select variable as LASSO or confirm prior information as Bayesian estimation.
Except the representation or model of neural network, it is the objective function that affects the usefulness of learning algorithms.

%For more on **loss function** see:
%
%* <https://blog.algorithmia.com/introduction-to-loss-functions/>;
%* <https://www.learnopencv.com/understanding-activation-functions-in-deep-learning/>;
%* <http://laid.delanover.com/activation-functions-in-deep-learning-sigmoid-relu-lrelu-prelu-rrelu-elu-softmax/>;
%* https://www.wikiwand.com/en/Robust_statistics
%* https://www.wikiwand.com/en/Huber_loss
%* <https://github.com/Kulbear/deep-learning-nano-foundation/wiki/ReLU-and-Softmax-Activation-Functions>;
%* <https://www.cs.umd.edu/~tomg/projects/landscapes/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Backpropagation, Training and Regularization}

\subsection{Backpropagation}

Automatic differentiation is the generic name for techniques that use the computational representation of a function to produce \textbf{analytic} values for the derivatives.
Automatic differentiation techniques are founded on the observation that 
any function, no matter how complicated, is evaluated by performing a sequence of simple elementary operations involving just one or two arguments at a time.
Backpropagation is one special case of automatic differentiation, i.e. \emph{reverse-mode automatic differentiation}.

The backpropagation procedure to compute the gradient of an objective function with respect to the weights of a multilayer stack of modules is 
nothing more than a practical application of the \textbf{chain rule for derivatives}.
The key insight is that the derivative (or gradient) of the objective with respect to the input of a module can be computed 
by working backwards from the gradient with respect to the output of that module (or the input of the subsequent module).
The backpropagation equation can be applied repeatedly to
propagate gradients through all modules, starting from the output at the top (where the network produces its prediction) all the way to the bottom (where the external input is fed).
Once these gradients have been computed, it is straightforward to compute the gradients with respect to the weights of each module.

%***
Suppose that $f(x)={\sigma}\circ(WH + b)$, where $H=\sigma\circ(W_4H_3 + b_4)$, $H_3=\sigma\circ(W_3H_2 + b_3)$, $H_2=\sigma\circ(W_2H_1 + b_2), $ $H_1=\sigma\circ(W_1x + b_1)$,
we want to compute the gradient $L(x_0,d_0)=\|f(x_0)-d_0\|^{2}_2$ with respect to all weights $W_1,W_2,W_3,W$:

$$\frac{\partial L(x_0,d_0)}{\partial W_n^i}=\frac{\partial L(x_0,d_0)}{\partial f(x_0)}\frac{\partial f(x_0)}{\partial W_n^i}$$
 where $\forall i\in\{1,2,\dots,l_n\}, \forall\,n\in\{1,2,3,4\}$
and it is fundamental to compute the gradient with respect to the last layer as below.
\begin{itemize}
	\item the gradient of loss function with respect to the prediction function:
		$$\frac{\partial L(x_0,d_0)}{\partial f(x_0)}=2[f(x_0)-d_0],$$

	\item the gradient of each unit in prediction function with respect to the weight in the last layer:
		$$
		\frac{\partial f^{j}(x_0)}{\partial W^j}=
		\frac{\partial \sigma(W^jH+b^j)}{\partial W^j}=
		{\sigma}^{\prime}(W^jH+b^j) H \,\,\forall j\in\{1,2,\dots,l\},
		$$

	\item the gradient of prediction function with respect to the last hidden state:
		$$
		\frac{\partial f^{j}(x_0)}{\partial H}  =
		\frac{\partial \sigma(W^jH+b^j)}{\partial H}  =
		{\sigma}^{\prime}(W^jH+b^j) W^j \,\,\forall j\in\{1,2,\dots,l\},
		$$
where $f^{j}(x_0)$, $W^{j}$, $b^j$ and $\sigma^{\prime}(z)$is the $j$th element of $f(x_0)$, 
the $j$-th row of matrix $W$, the $j$th element of vector $b$ and $\frac{\mathrm{d}\sigma(z)}{\mathrm{d} z}$, respectively.
\end{itemize}


\begin{figure}[!htbp]
\centering
\includegraphics[width=0.6\tectwidth]{MLP.jpg}
\caption{The Architecture of Feedforward Neural Networks}/label{fig:MLP4}
\end{figure}

\begin{note}
  Each connection ,the black line, is attached with a weight parameter.
\end{note}
%***

Recall the chain rule with more variables:
$$\frac{\partial f(m(x_0),n(x_0))}{\partial x_0}
	=  \frac{\partial f(m(x_0),n(x_0))}{\partial m(x_0)}\frac{\partial m(x_0)}{\partial x_0} + 
	  \frac{\partial f(m(x_0),n(x_0))}{\partial n(x_0)}\frac{\partial n(x_0)}{\partial x_0}.$$

Similarly , we can compute following gradients:
$$\frac{\partial H^j}{\partial W_4^j}  =\frac{\partial \sigma(W_4^j H_3+b_4^j)}{\partial W_4^j}  =\sigma^{\prime}(W_4^j H_3+b_4^j)H_3    \qquad\forall j\in\{1,2,\dots,l\};$$
$$\frac{\partial H^j}{\partial H_3}    =\frac{\partial \sigma(W_4^j H_3+b_4^j)}{\partial H_3}  =\sigma^{\prime}(W_4^j H_3+b_4^j)W_4^j  \qquad\forall j\in\{1,2,\dots,l_4\};$$
$$\frac{\partial H_3^j}{\partial W_3^j}=\frac{\partial \sigma(W_3^j H_2+b_3^j)}{\partial W_3^j}  =\sigma^{\prime}(W_3^j H_2+b_3^j)H_2    \qquad\forall j\in\{1,2,\dots,l_3\};$$
$$\frac{\partial H_3^j}{\partial H_2}  =\frac{\partial \sigma(W_3^j H_2+b_3^j)}{\partial H_2}  =\sigma^{\prime}(W_3^j H_2+b_3^j)W_3^j  \qquad\forall j\in\{1,2,\dots,l_3\};$$
$$\frac{\partial H_2^j}{\partial W_2^j}=\frac{\partial \sigma(W_2^j H_1+b_2^j)}{\partial W_2^j}  =\sigma^{\prime}(W_2^j H_1+b_2^j)H_1    \qquad\forall j\in\{1,2,\dots,l_2\};$$
$$\frac{\partial H_2^j}{\partial H_1}  =\frac{\partial \sigma(W_2^j H_1+b_2^j)}{\partial H_1}    =\sigma^{\prime}(W_2^j H_1+b_2^j)W_2^j  \qquad\forall j\in\{1,2,\dots,l_2\};$$
$$\frac{\partial H_1^j}{\partial W_1^j}=\frac{\partial \sigma(W_1^j x_0+b_1^j)}{\partial W_1^j}  =\sigma^{\prime}(W_1^j x_0+b_1^j)x_0    \qquad\forall j\in\{1,2,\dots,l_1\}.$$

%***

The multilayer perceptron $f(x)$ can be written in a chain form:
\begin{align*}
X              & \stackrel{\sigma}{\to} H_1   \stackrel{\sigma}{\to} H_2   \stackrel{\sigma} \to H_3    \stackrel{\sigma} \to H_4    \stackrel{\sigma} \to H   \stackrel{\sigma}\to f(x)  \\
\mathbb{R}^{p} & \rightarrow \mathbb{R}^{l_1} \rightarrow \mathbb{R}^{l_2} \rightarrow \mathbb{R}^{l_3} \rightarrow \mathbb{R}^{l_4} \rightarrow \mathbb{R}^{l}\rightarrow \mathbb{R}^{o} \\
X              & \rightarrow W_1 X \rightarrow W_2H_1 \rightarrow W_3H_2 \rightarrow W_4H_3 \rightarrow WH \rightarrow y       
\end{align*}
while the backpropagation to compute the gradient is in the reverse order:
\begin{equation}
\frac{\partial y}{\partial W}\to \frac{\partial y}{\partial H}\to \frac{\partial H}{\partial W_4}
\to \frac{\partial H}{\partial H_3}\to \frac{\partial H_3}{\partial W_3}
\to \frac{\partial H_3}{\partial H_2}\to \frac{\partial H_2}{\partial W_2}
\to \frac{\partial H_2}{\partial W_1}\to \frac{\partial H_1}{\partial W_1}.
\end{equation}

In general, the gradient of any weight can be computed by *backpropagation* algorithm.
The first step is to compute the gradient of loss function with respect to the output $f(x_0)=y\in\mathbb{R}^{o}$, i.e.
$\frac{\partial L(x_0, d_0)}{\partial f(x_0)}=2(f(x_0)-d_0)=2(\sigma\circ(WH+b)-d_0)$, 
of which the $i$th element is $2(y^{i}-d_0^i)=2(\sigma(W^{i}H+b^{i})-d_0^{i})\,\forall i\{1,2,\dots,o\}$.
Thus
$$\frac{\partial L(x_0, d_0)}{\partial W^{i}} = \frac{\partial L(x_0, d_0)}{\partial y^{i}}\frac{\partial y^{i}}{\partial W^{i}}=2(y^{i}-d_0^i)\sigma^{\prime}(W^iH+b^i)H.$$
Thus we can compute all the gradients of $W$ columns. Note that $H$ has been computed through forwards propagation in that layer.

And $H=\sigma\circ(W_4H_3+b_3)$, of which the $i$th element is $H^{i}=\sigma(W_4 H_3 +b_4)^{i}=\sigma(W_4^{i} H_3+b_4^{i})$.

And we can compute the gradient of columns of $W_4$:
\begin{align}
\frac{\partial L(x_0,y_0)}{\partial W_4^i}
& = \sum_{j=1}^{o}
[\frac{\partial L(x_0,y_0)}{\partial f^j (x_0)}]
[\frac{\partial f^j (x_0)}{\partial z}]
\frac{\partial z}{\partial W_4^i}
& =\sum_{j=1}^{o}
[\frac{\partial L(x_0,y_0)}{\partial y^j}]
[\frac{\partial y^j}{\partial H^i}]
\frac{\partial H^i}{\partial W_4^i} \\
& = \color{aqua}{
\sum_{j=1}^{o} 
[\frac{\partial L}{\partial y^j}]
[\frac{\partial\, y^j}{\partial (W^jH+b^j)}]
[\frac{\partial (W^jH+b^j)}{\partial H^{i}}]
\frac{\partial (H^{i})}{\partial W_4^i} } \\
&= \sum_{j=1}^{l}[\frac{\partial L}{\partial y^j}]\,[\sigma^{\prime}(W^j H+b^j)\,W^{j,i}][\,\sigma^{\prime}(W^i_4 H_3+b^i_4)H_3],
\end{align}
where $W^{j,i}$ is the $i$th element of $j$th column in matrix $W$.

\begin{align}
\frac{\partial L(x_0,y_0)}{\partial W_3^i}
&=\sum_{j=1}^{o}
\frac{\partial L(x_0,y_0)}{\partial f^j (x_0)}
[\frac{\partial f^j (x_0)}{\partial z}]
\frac{\partial z}{\partial W_3^i} \\
&=\sum_{j=1}^{o}
\frac{\partial L}{\partial y^j }
[\frac{\partial y^j}{\partial H_3^i}]
\frac{\partial H_3^i}{\partial W_3^i}\\
&=\sum_{j=1}^{o}
\frac{\partial L}{\partial y^j }
[\sum_{k=1}\frac{\partial y^j}{\partial H^k} \frac{\partial H^k}{\partial H_3^i}]
\frac{\partial H_3^i}{\partial W_3^i}
\end{align}
where all the partial derivatives or gradients have been computed or accessible. 
It is nothing except to add or multiply these values in the order when we compute the weights of hidden layer.


\begin{align}
\frac{\partial L(x_0,y_0)}{\partial W_2^i}
& =\sum_{j=1}^{o}
\frac{\partial L(x_0,y_0)}{\partial f^j (x_0)}
[\frac{\partial f^j (x_0)}{\partial z}]
\frac{\partial z}{\partial W_2^i}      \\
& =\sum_{j=1}^{l}
\frac{\partial L}{\partial y^j }
[\frac{\partial y^j}{\partial H_2^i}]
\frac{\partial H_2^i}{\partial W_2^i}  \\
& =\sum_{j=1}^{o}
\frac{\partial L}{\partial y^j }
\{\sum_{k=1}\frac{\partial y^j}{\partial H^k}
[\sum_{m}\frac{\partial H^k}{\partial H_3^m}\frac{\partial H_3^m}{\partial H_2^i}]\}
\frac{\partial H_2^i}{\partial W_2^i}
\end{align}

And the gradient of the first layer is computed by
\begin{align}
\frac{\partial L(x_0,y_0)}{\partial W_1^i}
& =\sum_{j}\frac{\partial L(x_0,y_0)}{\partial y^j}\frac{\partial y^j}{\partial z}\frac{\partial z}{\partial W_1^i}                              \\
& =\sum_{j}\frac{\partial L}{\partial y^j}
[\sum_{k}\frac{\partial y^j}{\partial H^k}
\sum_{m}\frac{\partial H^k}{\partial H_3^m}
\sum_{n}\frac{\partial H_3^k}{\partial H_2^n}
\sum_{r}\frac{\partial H_2^n}{\partial H_1^r}]
\frac{\partial H_1^i}{\partial W_1^i}.
\end{align}
See more information on backpropagation in the following list
\begin{itemize}
\item \href{http://www.offconvex.org/2016/12/20/backprop/}{Back-propagation, an introduction at offconvex.org};
%* [Backpropagation on Wikipedia](https://www.wikiwand.com/en/Backpropagation);
\item \href{https://www.wikiwand.com/en/Automatic_differentiation}{Automatic differentiation on Wikipedia};
\item \href{https://brilliant.org/wiki/backpropagation/}{backpropagation on brilliant};
%* An introduction to automatic differentiation at <https://alexey.radul.name/ideas/2013/introduction-to-automatic-differentiation/>;
%* Reverse-mode automatic differentiation: a tutorial at <https://rufflewind.com/2016-12-30/reverse-mode-automatic-differentiation>.
\item \href{http://people.idsia.ch/~juergen/who-invented-backpropagation.html}{Who invented backpropagation ?};
%* [Autodiff Workshop](https://autodiff-workshop.github.io/)
%* [如何直观地解释 backpropagation 算法？ - 景略集智的回答 - 知乎](https://www.zhihu.com/question/27239198/answer/537357910)
%* The chapter 2 *How the backpropagation algorithm works* at the online book <http://neuralnetworksanddeeplearning.com/chap2.html>
\item For more information on automatic differentiation see the book \emph{Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Second Edition} 
      by Andreas Griewank and Andrea Walther_ at \url{https://epubs.siam.org/doi/book/10.1137/1.9780898717761}.
\end{itemize}
%![](http://ai.stanford.edu/~tengyuma/forblog/weight5.jpg)

\subsection{Training Methods}

The training is to find the optimal parameters of the model based on the \textbf{training data set}. 
The training methods are usually based on the gradient of cost function as well as back-propagation algorithm  in deep learning.
See \emph{Stochastic Gradient Descent} in \emph{Numerical Optimization} for details.
In this section, we will talk other optimization tricks such as \emph{Normalization}.
\ebgin{itemize}
  \item \href{https://arxiv.org/abs/1803.08494}{Group Normalization} at \url{https://arxiv.org/abs/1803.08494}.
  \item \href{https://blog.paperspace.com/busting-the-myths-about-batch-normalization/}{Busting the myth about batch normalization at paperspace.com}.
  \item The original paper \emph{Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift} at \url{https://arxiv.org/pdf/1502.03167.pdf}.
\end{itemize}


\subsection{Regularization}

In mathematics, statistics, and computer science, particularly in the fields of machine learning and inverse problems, 
regularization is a process of introducing additional information 
in order to solve an ill-posed problem or to prevent over-fitting.
In general, regularization is a technique that applies to objective functions in ill-posed optimization problems.
It changes the objective function or more generally the optimization procedure. 
However, it is not crystal clear that what is the relationship between the optimization techniques and generalization ability.
See the following links for more information on optimization and generalization.

\begin{itemize}
  \item \url{https://www.inference.vc/sharp-vs-flat-minima-are-still-a-mystery-to-me/}
  \item \url{https://arxiv.org/abs/1703.04933}
  \item \url{https://arxiv.org/abs/1810.05369}
  \item \url{https://blog.csdn.net/xzy_thu/article/details/80732220}
  \item \url{http://www.offconvex.org/2017/12/08/generalization1/}
  \item \url{http://www.offconvex.org/2018/02/17/generalization2/}
%* http://www.offconvex.org/2017/03/30/GANs2/
%* https://machinelearningmastery.com/blog/
  \item \url{http://lcsl.mit.edu/courses/regml/regml2016/}
  \item \url{https://chunml.github.io/ChunML.github.io/tutorial/Regularization/}
  \item \url{http://www.mit.edu/~9.520/fall16/}
  \item \url{https://arxiv.org/pdf/1506.02142.pdf}
\end{itemize}

\subsubsection{Parameter norm penalty}

The $\ell_2$ norm penalty  is to add the squares of $\ell_2$ norm of parameters to the objective function $J(\theta)$ to reduce the parameters(or weights) 
as shown in ridge regression with regular term coefficient $\lambda$, i.e.
$ J(\theta)+\lambda {\|\theta\|}_{2}^{2}.$
Suppose  $E(\theta)=J(\theta)+\lambda {\|\theta\|}_{2}^{2}$, the gradient descent take approximate (maybe inappropriate)  form
$$ \theta=\theta-\eta\frac{\partial E(\theta)}{\partial \theta}=\theta -\eta\frac{\partial J(\theta)}{\partial \theta}-2\eta\lambda \theta$$
thus
$$ \frac{\partial J(\theta)}{\partial \theta}=-2\lambda\theta\implies J(\theta)=e^{-2\lambda \theta}. $$

If we want to find the minima of $E(\theta)$, $\theta$ will decay to $0$.
It extends to the following iterative formula:
$$\theta^{t+1} = (1-\lambda)\theta^{t}-\alpha_{t}\frac{\partial J(\theta^{t})}{\partial \theta},$$
where $\lambda$  determines how you trade off the original cost $J(\theta)$ with the large weights penalization.
The new term $\lambda$ coming from the regularization causes the weight to decay in proportion to its size.

\begin{itemize}
  \item \url{https://stats.stackexchange.com/questions/70101/neural-networks-weight-change-momentum-and-weight-decay};
  \item \url{https://metacademy.org/graphs/concepts/weight_decay_neural_networks};
  \item \url{https://machinelearningmastery.com/how-to-reduce-overfitting-in-deep-neural-networks-with-weight-constraints-in-keras/};
  \item \url{https://machinelearningmastery.com/how-to-reduce-overfitting-in-deep-learning-with-weight-regularization/}.
\end{itemize}

The $\ell_1$ norm penalty is also used in deep learning as in **LASSO**. It is to solve the following optimization problem:
  $$\min_{\theta}J(\theta)+\lambda{\|\theta\|}_1, $$
where $\lambda$ is a hyperparameter. Sparsity  brings to the model as shown as in **LASSO**.

\subsubsection{Early stop}

Its essential is to make a balance in memorization and generalization.
Early stopping is to stop the procedure before finding the minima of cost in training data. 
It is one direct application of \emph{cross validation}.

\begin{itemize}
	\item \url{ https://www.wikiwand.com/en/Early_stopping};
	\item \url{https://www.wikiwand.com/en/Cross-validation_(statistics)};
	\item \url{https://machinelearningmastery.com/how-to-stop-training-deep-neural-networks-at-the-right-time-using-early-stopping/};
	\item \url{https://machinelearningmastery.com/early-stopping-to-avoid-overtraining-neural-network-models/}.
\end{itemize}


\subsubsection{Dropout}

It is to cripple the connections stochastically, which  is often used in visual tasks. 
See the original paper \href{http://www.cs.toronto.edu/~hinton/absps/JMLRdropout.pdf}{Dropout: A Simple Way to Prevent Neural Networks from Overfitting}.

\iffalse 
* https://www.zhihu.com/question/24529483
* https://www.jiqizhixin.com/articles/2018-11-10-7
* https://www.jiqizhixin.com/articles/112501
* https://www.jiqizhixin.com/articles/2018-08-27-12
* https://yq.aliyun.com/articles/68901
* https://www.wikiwand.com/en/Regularization_(mathematics)
* [CNN tricks](http://lamda.nju.edu.cn/weixs/project/CNNTricks/CNNTricks.html)
* https://www.jeremyjordan.me/deep-neural-networks-preventing-overfitting/
* https://www.doc.ic.ac.uk/~nd/surprise_96/journal/vol4/cs11/report.html#An%20engineering%20approach
* https://machinelearningmastery.com/dropout-for-regularizing-deep-neural-networks/
\fi

\subsubsection{Data Augmentation}

Data augmentation is to augment the training datum specially in visual recognition.
\textbf{Overfitting} in supervised learning is data-dependent. In other words, the model may generalize better if the data set is more diverse.
It is to collect more datum in the statistical perspective.

%* [The Effectiveness of Data Augmentation in Image Classification using Deep Learning](http://cs231n.stanford.edu/reports/2017/pdfs/300.pdf)
%* http://www.cnblogs.com/love6tao/p/5841648.html

%***

%|Feed forward and Propagate backwards|
%|:----------------------------------:|
%|![TF](http://beamandrew.github.io//images/deep_learning_101/tensors_flowing.gif)|

\chapter{ElegantBook 设置说明}


\section{编译方式}

本模板基于基础的 book 文类，所以 book 的选项对于本模板也是有效的。默认编码为 UTF-8，推荐使用 \TeX{} Live 编译。作者编写环境为 Win10(64bit) + \TeX{} Live 2018。由于使用的是 \texttt{ctex} 宏包，所以支持 \texttt{pdflatex} 以及 \texttt{xelatex} 编译。


\section{选项设置}
本文特殊选项设置共有 2 类，分为 {\color{main}主题颜色}设置 以及 {\color{main}章标题显示风格}设置。

第 1 类为{\color{main}主题颜色}设置，内置 3 组颜色主题，分别为 \verb|green|（默认），\verb|cyan|，\verb|blue|，另外还有一个自定义的选项 \verb|nocolor|，用户\textbf{必须}在使用模板的时候选择某个颜色主题或选择 \verb|nocolor| 选项。调用颜色主题 \verb|green| 的方法为 \verb|\documentclass[green]{elegantbook}|或者使用 \verb|\documentclass[color=green]{elegantbook}|。需要改变颜色的话请选择 \verb|nocolor| 选项或者使用 \verb|color=none|，然后在导言区定义 main、second、third 颜色，具体的方法如下：
\begin{verbatim}
\definecolor{main}{RGB}{70,70,70}    %定义 main 颜色值
\definecolor{second}{RGB}{115,45,2}    %定义 second 颜色值
\definecolor{third}{RGB}{0,80,80}     %定义 third 颜色值
\end{verbatim}

\begin{table}[htp]
\centering
\begin{tabular}{ccccc}
\toprule
	  & green & cyan & blue & 主要使用的环境\\
\midrule
main & \makecell{{\color{main1}\rule{1cm}{1cm}}}& \makecell{{\color{main2}\rule{1cm}{1cm}}}&\makecell{ {\color{main3}\rule{1cm}{1cm}}}& definition \\

second &\makecell{ {\color{second1}\rule{1cm}{1cm}}}& \makecell{{\color{second2}\rule{1cm}{1cm}}}&\makecell{ {\color{second3}\rule{1cm}{1cm}}}&theorem \ lemma \ corollary\\

third &\makecell{ {\color{third1}\rule{1cm}{1cm}}}& \makecell{{\color{third2}\rule{1cm}{1cm}}}&\makecell{ {\color{third3}\rule{1cm}{1cm}}}&proposition\\
\bottomrule
\end{tabular}
\caption{ElegantBook 模板中的三套颜色主题\label{tab:color thm}}
\end{table}

第 2 类为{\color{main} 章标题显示风格}，包含 \verb|hang|（默认）与 \verb|display| 两种风格，区别在于章标题单行显示（\verb|hang|）与双行显示（\verb|display|），本说明使用了 \verb|hang|。调用方式为 \verb|\documentclass[hang]{elegantbook}| 或者 \verb|\documentclass[titlestyle=hang]{elegantbook}|。

综合起来，同时调用三个选项使用 \verb|\documentclass[color=X,titlestyle=Y]{elegantbook}|。其中 \verb|X| 可以选择 \verb|green|,\verb|cyan|,\verb|blue|,\verb|none|；\verb|Z| 可以选择 \verb|hang| 或者 \verb|display|。

\section{数学环境简介}

在我们这个模板中，定义了三大类环境

\begin{enumerate}[noitemsep]
\item 定理类环境，包含标题和内容两部分。根据格式的不同分为3种
   \begin{itemize}[noitemsep]
      \item {\color{main}\bfseries definition} 环境，含有一个可选项，编号以章节为单位，颜色为 {\color{main}main}；
      \item {\color{second}\bfseries theorem、lemma、corollary} 环境，颜色为主颜色 {\color{second}second}，编号均以章节为单位；
      \item {\color{third}\bfseries proposition} 环境，含有一个可选项，编号以章节为单位，颜色为 {\color{third}{third}}。
   \end{itemize}
\item 示例类环境，有 \textbf{example、exercise} 环境，自动编号，编号以章节为单位。
\item 证明类环境，有 \textbf{proof、note} 环境，特点是，有引导符或者结尾符，\textbf{note} 环境有引导符号，\textbf{proof} 环境有证明完毕标志。
\item 结论类环境，有\textbf{conclusion、assumption、property，remark、solution} 环境，三者均以粗体的引导词为开头，和普通段落格式一致。
\end{enumerate}

\section{可编辑的字段}
在模板中，可以编辑的字段分别为作者 \verb|\author|、邮箱 \verb|\email|、中文标题 \verb|\zhtitle|、中文标题结尾 \verb|\zhend|、英文标题\verb| \entitle|、英文标题结尾 \verb|\enend|、名言 \verb|\myquote|、版本号 \verb|\version|。并且，可以根据自己的喜好把封面水印效果的 \verb|cover.pdf| 替换掉，以及封面中用到的 \verb|logo.pdf|。

\chapter{ElegantBook 写作示例}

\section{Economics and Differentiable Function}

Economics focuses on the behaviour and interactions of economic agents and how economies work. Microeconomics analyzes basic elements in the economy, including individual agents and markets, their interactions, and the outcomes of interactions. Individual agents may include, for example, households, firms, buyers, and sellers. Macroeconomics analyzes the entire economy (meaning aggregated production, consumption, savings, and investment) and issues affecting it, including unemployment of resources (labour, capital, and land), inflation, economic growth, and the public policies that address these issues (monetary, fiscal, and other policies). See glossary of economics.
\begin{align*}
   &\max(\min)\quad \mathbb{E}\int_{t_0}^{t_1}f(t,x,u)\,dt\\
   &\quad\mbox{s.t.} \quad dx=g(t,x,u)dt+\sigma(t,x,u)dz\\
   &\quad \hspace{2.em} k(0)=k_0\;\text{given}
\end{align*}

where $z$ is stochastic process or white noise or wiener process.

Other broad distinctions within economics include those between positive economics, describing "what is", and normative economics, advocating "what ought to be"; between economic theory and applied economics; between rational and behavioural economics; and between mainstream economics and heterodox economics.

Economic analysis can be applied throughout society, in business, finance, health care, and government. Economic analysis is sometimes also applied to such diverse subjects as crime, education, the family, law, politics, religion, social institutions, war, science, and the environment.

\begin{definition}{Differenzierbarkeit}{diff}
Eine Funktion $f:~I\to\mathbb{R}$ auf einem Intervall $I$ hei\ss{}t in $x_0\in I$ differenzierbar oder linear approximierbar, wenn der Grenzwert
   \begin{equation*}
      \lim\limits_{x\to x_0}\frac{f(x)-f(x_0)}{x-x_0}=
      \lim\limits_{h\to 0}\frac{f(x_0+h)-f(x_0)}{h}
   \end{equation*}
existiert. Bei Existenz hei\ss{}t dieser Grenzwert Ableitung oder Differentialquotient von $f$ in $x_0$ und man schreibt f\"{u}r ihn
   \begin{equation*}
      f'(x_0)\quad\text{oder}\quad\frac{df}{dx}(x_0).
   \end{equation*}
\end{definition}

The discipline was renamed in the late 19th century, primarily due to Alfred Marshall, from "political economy" to "economics" as a shorter term for "economic science". At that time, it became more open to rigorous thinking and made increased use of mathematics, which helped support efforts to have it accepted as a science and as a separate discipline outside of political science and other social sciences.

\begin{example}
$E$ and $F$ be two events such that $P(E)=
P(F)=1/2$, and $P(E\cap F)=1/3$, let $\mathscr{F}=\sigma(Y)$,  $X$ and $Y$ be the indicate function of $E$ and $F$ respectively. How to compute $\mathbb{E}[ X\mid \mathscr{F} ]$?
\end{example}

Some subsequent comments criticized the definition as overly broad in failing to limit its subject matter to analysis of markets. From the 1960s, however, such comments abated as the economic theory of maximizing behaviour and rational-choice modelling expanded the domain of the subject to areas previously treated in other fields. There are other criticisms as well, such as in scarcity not accounting for the macroeconomics of high unemployment.

\begin{exercise}
let $S=l^\infty=\big\{(x_n)\mid \exists\, M \text{ such that } \forall n, |x_n|\leq M,x_n\in \mathbb{R}\big\}$, $\rho_{\infty}(x,y)=\sup\limits_{n\geq 1}|x_n-y_n|$, show that $\big(l^\infty,\rho_{\infty}\big)$ is complete.
\end{exercise}

\begin{theorem}{Mittelwertsatz f\"{u}r $n$ Variable}{31}
Es sei $n\in\mathbb{N}$, $D\subseteq\mathbb{R}^n$ eine offene Menge und $f\in C^{1}(D,\mathbb{R})$. Dann gibt es auf jeder Strecke $[x_0,x]\subset D$ einen Punkt $\xi\in[x_0,x]$, so dass gilt
\begin{equation*}
   f(x)-f(x_0) = \operatorname{grad} f(\xi)^{\top}(x-x_0)
\end{equation*}
\end{theorem}

\begin{note}
在本模板中，引理（lemma），推论（corollary ）的样式和定理的样式一致，包括颜色，仅仅只有计数器的设置不一样。在这个例稿中，我们将不给出引理推论的例子。
\end{note}


Gary Becker, a contributor to the expansion of economics into new areas, describes the approach he favours as "combin[ing the] assumptions of maximizing behaviour, stable preferences, and market equilibrium, used relentlessly and unflinchingly."One commentary characterizes the remark as making economics an approach rather than a subject matter but with great specificity as to the "choice process and the type of social interaction that [such] analysis involves."


Economic efficiency measures how well a system generates desired output with a given set of inputs and available technology. Efficiency is improved if more output is generated without changing inputs, or in other words, the amount of "waste" is reduced. A widely accepted general standard is Pareto efficiency, which is reached when no further change can make someone better off without making someone else worse off.

\begin{proposition}{最优性原理}{max}
如果 $u^*$ 在 $[s,T]$ 上为最优解，则 $u^*$在$[s,T]$ 任意子区间都是最优解，假设区间为 $[t_0,t_1]$ 的最优解为 $u^*$ ，则 $u(t_0)=u^{*}(t_0)$，即初始条件必须还是在 $u^*$ 上。
\end{proposition}

Microeconomics examines how entities, forming a market structure, interact within a market to create a market system. These entities include private and public players with various classifications, typically operating under scarcity of tradable units and light government regulation. The item traded may be a tangible product such as apples or a service such as repair services, legal counsel, or entertainment.

\begin{corollary}{}{}
假设 $V(\cdot,\cdot)$ 为值函数，则根据最大值原理~\ref{thm:31}，有如下推论
\[
V(k,z)=\max\Big\{u\big(zf(k)-y\big)+\beta \mathbb{E}V(y,z^\prime)\Big\}
\]
\end{corollary}

\begin{proof}
因为 $y^*=\alpha\beta z k^\alpha$，$V(k,z)=\alpha/1-\alpha\beta\ln k_0+1/1-\alpha\beta \ln z_0+\Delta$。
   \begin{align*}
   \text{右边}&=\Big\{u\big(zf(k)-y\big)+\beta \mathbb{E}V(y,z^\prime)\Big\}\\
   &=\ln(zk^\alpha-\alpha\beta zk^\alpha)+\beta\mathbb{E}\Big[\frac{\alpha}{1-\alpha\beta}\ln y+\frac{1}{1-\alpha\beta}\ln z^\prime+\Delta\Big]\\
   &=\ln(1-\alpha\beta)zk^\alpha+\beta\Big\{\mathbb{E}\big[\frac{\alpha}{1-\alpha\beta}\ln \alpha\beta z k^\alpha\big]+\frac{1}{1-\alpha\beta}\mathbb{E}[\ln z^\prime]+\Delta\Big\}
   \end{align*}
利用 $\mathbb{E}[\ln z^\prime]=0$，并将对数展开得
   \begin{align*}
   \text{右边}&=\ln (1-\alpha\beta)+\ln z+\alpha\ln k+\frac{\alpha\beta}{1-\alpha\beta}\big[\ln \alpha\beta+\ln z+\alpha\ln k\big]+\frac{\beta}{1-\alpha\beta}\mu+\beta \Delta\\
   &=\frac{\alpha}{1-\alpha\beta}\ln k+\frac{1}{1-\alpha\beta}\ln z+\Delta
   \end{align*}
所以 $\text{左边}=\text{右边}$，证毕。
\end{proof}



\begin{property}
Properties of Cauchy Sequence
\begin{enumerate}[noitemsep]
\item $\{x_k\}$ is cauchy sequence then $\{x_k^i\}$ is cauchy sequence.
\item $x_k\in \mathbb{R}^n$, $\rho(x,y)$ is Euclidean, then cauchy is equivalent to convergent, $(\mathbb{R}^n,\rho)$ metric space is complete.
\end{enumerate}
\end{property}

\begin{note}
conclusion、assumption、property，remark、solution 的环境效果是一样的。
\end{note}

Various market structures exist. In perfectly competitive markets, no participants are large enough to have the market power to set the price of a homogeneous product. In other words, every participant is a "price taker" as no participant influences the price of a product. In the real world, markets often experience imperfect competition.

Scarcity is represented in the figure by people being willing but unable in the aggregate to consume beyond the PPF (such as at $X$) and by the negative slope of the curve.[32] If production of one good increases along the curve, production of the other good decreases, an inverse relationship. This is because increasing output of one good requires transferring inputs to it from production of the other good, decreasing the latter.

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=0.6\textwidth]{mpg.png}
	\caption{The Relationship between MPG and Weight\label{fig:mpg}}
\end{figure}

\begin{definition}{Contraction mapping}{mapping}
$(S,\rho)$ is the metric space, $T: S\to S$, If there exists $\alpha\in(0,1)$ such that for any $x$ and $y\in S$, the distance
\begin{equation}
\rho(Tx,Ty)\leq \alpha\rho(x,y)
\end{equation}
Then $T$ is a {\color{main} contraction mapping}.
\end{definition}

\begin{remark}
\begin{enumerate}[noitemsep]
\item $T:S\to S$, where $S$ is a metric space, if  for any $x,y\in S$, $\rho(Tx,Ty)<\rho(x,y)$ is not contraction mapping.
\item Contraction mapping is continuous map.
\end{enumerate}
\end{remark}


\begin{conclusion}
In theory, in a free market the aggregates (sum of) of quantity demanded by buyers and quantity supplied by sellers may reach economic equilibrium over time in reaction to price changes; in practice, various issues may prevent equilibrium, and any equilibrium reached may not necessarily be morally equitable. For example, if the supply of healthcare services is limited by external factors, the equilibrium price may be unaffordable for many who desire it but cannot pay for it.
\end{conclusion}

\section{Bibliography}
This template uses Bib\TeX{} to generate the bibliography, the default bibliography style is \verb|aer|. ~\cite{Chen2018} use data from a major peer-to-peer lending marketplace in China to study whether female and male investors evaluate loan performance differently. You can add bib items (from Google Scholar, Mendeley, EndNote, and etc.) to \verb|reference.bib| file, and cite the bibkey in the \verb|tex| file.


\nocite{EINAV2010,Havrylchyk2018}

\bibliographystyle{aer}
\bibliography{reference}

\addcontentsline{toc}{chapter}{参考文献}


\end{document}