Merge branch 'master' of https://github.com/CompPhysics/MachineLearning
This commit is contained in:
Christian Forssen
@@ -119,6 +119,249 @@ process then follows the same structure as the VMC method.
|
||||
|
||||
|
||||
|
||||
A common thread in many unsupervised learning tasks is accurately
|
||||
representing the underlying probability distribution from which a
|
||||
dataset is drawn. Unsupervised learning of high-dimensional, complex
|
||||
distributions presents a new set of technical and computational
|
||||
challenges that are different from those we encountered in a
|
||||
supervised learning setting. When dealing with complicated probability
|
||||
distributions, it is often much easier to learn the \emph{relative
|
||||
weights} of different states or data points (ratio of probabilities),
|
||||
than \emph{absolute} probabilities. In physics, this is the familiar
|
||||
statement that the weights of a Boltzmann distribution are much easier
|
||||
to calculate than the partition function. The relative probability of
|
||||
two configurations, $\bd{x}_1$ and $\bd{x}_2$, are proportional to the
|
||||
difference between their Boltzmann weights \be {p(\mathbf{x}_1) \over
|
||||
p(\mathbf{x}_2)} =e^{-\beta
|
||||
\left(E(\mathbf{x}_1)-E(\mathbf{x}_2)\right)}, \ee where as is usual
|
||||
in statistical mechanics $\beta$ is the inverse temperature and
|
||||
$E(\mathbf{x}; \theta)$ is the energy of state $\mathbf{x}$ given some
|
||||
parameters (couplings) $\theta$ . However, calculating the absolute
|
||||
weight of a configuration requires knowledge of the partition function
|
||||
\be Z_p =\mathrm{Tr}_{\mathbf{x} }e^{-\beta E(\mathbf{x})}, \ee (where
|
||||
the trace is taken over all possible configurations $\mathbf{x}$)
|
||||
since \be p(\mathbf{x}) ={ e^{-\beta E(\mathbf{x}_1)} \over Z_p}. \ee
|
||||
In general, calculating the partition function $Z_p$ is analytically
|
||||
and computationally intractable.
|
||||
|
||||
For example, for the Ising model with $N$ binary spins, the trace involves calculating a sum over $2^N$ terms, which is a difficult task for most energy functions. For this reason, physicists (and machine learning scientists) have developed various numerical and computational methods for evaluating such partition functions. One approach is to use Monte-Carlo based methods to draw samples from the underlying distribution (this can be done knowing only the relative probabilities) and then use these samples to numerically estimate the partition function. This is the philosophy behind powerful methods such as Markov Chain Monte Carlo (MCMC) \cite{andrieu2003introduction} and annealed importance sampling \cite{neal1998view} which are widely used in both the statistical physics and machine learning communities. An alternative approach -- which we focus on here -- is to approximate the the probability distribution $p(\mathbf{x})$ and partition function using a ``variational distribution'' $q(\mathbf{x}; \theta_q)$ whose partition function we can calculate exactly. The variational parameters $\theta_q$ are chosen to make the variational distribution as close to the true distribution as possible (how this is done is the focus of much of this section).
|
||||
|
||||
One of the most-widely applied examples of a variational method in statistical physics is Mean-Field Theory (MFT). MFT can be naturally understood as a procedure for approximating the true distribution of the system by a factorized distribution. The deep connection between MFT and variational methods is discussed below. These variational MFT methods have been extended to understand more complicated spin models (also called graphical models in the ML literature) and form the basis of powerful set of techniques that go under the name of Belief Propagation and Survey Propagation \cite{yedidia2003understanding, mackay2003information, wainwright2008graphical}.
|
||||
|
||||
Variational methods are also widely used in ML to approximate complex probabilistic models. For example, below we show how the Expectation Maximization (EM) procedure, which we discussed in the context of Gaussian Mixture Models for clustering, is actually a general method that can be derived for any latent (hidden) variable model using a variational procedure \cite{neal1998view}. This section serves as an introduction to this powerful class of variational techniques. For readers interested in an in-depth discussion on variational inference for probabilistic graphical models, we recommend the great treatise written by Michael I. Jordan and others\cite{jordan1999introduction}, the more physics oriented discussion in \cite{yedidia2001idiosyncratic, yedidia2003understanding}, as well as David MacKay's outstanding book \cite{mackay2003information}.
|
||||
|
||||
%
|
||||
%In applying variational methods, one often approximate the complex model of interest with a proxy that depends on some ``variational parameters". One then optimize an objective function with respect to these parameters to improve the accuracy of approximation. For example, in both physics and machine learning we usually have to deal with probability distributions of the form
|
||||
%\begin{equation}
|
||||
%p(\bd{s}|\beta,\bd{J}) = \frac{1}{Z(\bd{J}|\beta)}\mathrm e^{-\beta E(\bd{s},\bd{J})}
|
||||
%\end{equation}
|
||||
%where $E(\bd{s},\bd{J})$ is the ``energy'' of a configuration $\bd{s}$, and the normalization constant, or partition function, is given by $Z(\bd{J}|\beta) = \sum_{\bd{s}}\mathrm e^{-\beta E(\bd{s},\bd{J})}$. In this case, $\bd{J}$ is interpreted as the parameters while $\bd{s}$ as the data.
|
||||
%
|
||||
%It is not difficult to see that in many cases the partition function could be difficult to compute. This is problematic because in the physics context, it allows us to derive all thermodynamic quantities of interest, and in the machine learning context, normalization is important for inference on, say, posterior distribution . As we stressed, one can tackle this problem with approaches based on random sampling such as MCMC or more refined algorithm such as annealed importance sampling\cite{neal1998view}, here we'll rather be concerned with variational methods that bear a close resemblance to the variational free energy approach of Feynmann and Bogoliubov.
|
||||
%Therefore, the goal is to find approximate probability distributions $q(\bd{s},\bd{\theta})$ which attempts to capture the relevant features of the true distribution $p(\bd{s}|\beta,\bd{J})$ and, at the same time, are simple enough to allow the calculation of the corresponding partition function. To this end, we first recall variational mean-field theory using the Ising model as an illustrative example.
|
||||
|
||||
\subsection{\label{subsec:MFT_Ising}Variational mean-field theory for the Ising model}
|
||||
|
||||
Ising models are a major paradigm in statistical physics. Historically introduced to study magnetism, it was quickly realized that their predictive power applies to a variety of interacting many-particle systems. Ising models are now understood to serve as minimal models for complex phenomena such as certain classes of phase transitions. In the Ising model, degrees of freedom called spins assume discrete, binary values, e.g.~$s_i=\pm 1$. Each spin variable $s_i$ lives on a lattice (or, in general, a graph), the sites of which are labeled by $i=1,2\dots, N$. Despite the extreme simplicity relative to real-world systems, Ising models exhibit a high level of intrinsic complexity, and the degrees of freedom can become correlated in sophisticated ways. Often, spins interact spatially locally, and respond to externally applied magnetic fields.
|
||||
|
||||
A spin configuration $\bd{s}$ specifies the values $s_i$ of the spins at every lattice site. We can assign an ``energy'' to every such configuration
|
||||
\begin{equation}
|
||||
E(\bd{s},\bd{J}) = -\frac{1}{2}\sum_{i,j}J_{ij}s_i s_j - \sum_i h_i s_i,
|
||||
\end{equation}
|
||||
where $h_i$ is a local magnetic field acting on the spin $s_i$, and $J_{ij}$ is the interaction strength between the spins $s_i$ and $s_j$. In textbook examples, the coupling parameters $\bd{J}=(J,h)$ are typically uniform or, in studies of disordered systems, drawn from some probability distribution (i.e. quenched disorder).
|
||||
|
||||
The probability of finding the system in a given spin configuration at temperature $\beta^{-1}$ is given by
|
||||
\begin{eqnarray}
|
||||
p(\bd{s}|\bd{J}) &=& \frac{1}{Z_p(\bd{J})}\mathrm e^{-\beta E(\bd{s},\bd{J})},\nonumber\\
|
||||
Z_p(\bd{J})&=&\sum_{\{s_i=\pm 1\}} \mathrm e^{-\beta E(\bd{s},\bd{J})},
|
||||
\end{eqnarray}
|
||||
with $\sum_{\{s_i=\pm 1\}}$ denoting the sum over all possible configurations of the spin variables. We write $Z_p$ to emphasize that this is the partition function corresponding to the probability distribution $p(\bd{s})$, which will become important later. For a fixed number of lattice sites $N$, there are $2^N$ possible configurations, a number that grows exponentially with the system size. Therefore, it is not in general feasible to evaluate the partition function $Z_p(\bd{J})$ in closed form. This represents a major practical obstacle for extracting predictions from physical theories since the partition function is directly related to the free-energy through the expression
|
||||
\begin{eqnarray}
|
||||
\label{eq:exact_F}
|
||||
\beta F_p(\bd{J})&=&-\log Z_p(\bd{J})= \beta \langle E(\bd{s},\bd{J})\rangle_p - H_p,
|
||||
\end{eqnarray}
|
||||
with
|
||||
\begin{equation}
|
||||
\label{eq:entropy}
|
||||
H_p = -\sum_{\{s_i=\pm 1\}}p(\bd{s}|\bd{J})\log p(\bd{s}|\bd{J})
|
||||
\end{equation}
|
||||
the entropy of the probability distribution $p(\bd{s}|\bd{J})$.
|
||||
|
||||
|
||||
Even though the true probability distribution $p(\bd{s}|\beta,\bd{J})$ may be a very complicated object, we can still make progress by approximating $p(\bd{s}|\beta,\bd{J})$ by a \emph{variational distribution} $q(\bd{s},\bd{\theta})$ which captures the essential features of interest, with $\theta$ some parameters that define our variational ansatz. The name variational distribution comes from the fact that we are going to vary the parameters $\bd{\theta}$ to make $q(\bd{s},\bd{\theta})$ as close to $p(\bd{s}|\beta,\bd{J})$ as possible. The functional form of $q(\bd{s},\bd{\theta})$ is based on an ``educated guess'', which oftentimes comes from our intuition about the problem. We can also define a variational free-energy
|
||||
\be
|
||||
\beta F_q(\bd{\theta}, \bd{J}) = \beta \langle E(\bd{s},\bd{J})\rangle_q - H_q,
|
||||
\label{eq:main_var_eq}
|
||||
\ee
|
||||
where $ \langle E(\bd{s},\bd{J})\rangle_q$ is the expectation value of the energy corresponding to the distribution $p(\bd{s})$ with respect to the distribution $q(\bd{s},\bd{\theta})$, and $H_q$ is the entropy of $q(\bd{s},\bd{\theta})$.
|
||||
|
||||
Before proceeding further, it is helpful to introduce a new quantity: the Kullback-Leibler divergence (KL-divergence or relative entropy) between two distributions $p(\mathbf{x})$ and $q(\mathbf{x})$. The KL-divergence measures the dissimilarity between the two distributions and is given by
|
||||
\be
|
||||
D_{KL}(q\| p)=\mathrm{Tr}_{\mathbf{x}} q(\mathbf{x})\log{q(\mathbf{x}) \over p(\mathbf{x})},
|
||||
\ee
|
||||
which is the expectation w.r.t.~$q$ of the logarithmic difference between the two distributions $p$ and $q$. Two important properties of the KL-divergence are (i) positivity: $D_{KL}(p\| q)\geq 0$ with equality if and only if $p=q$ (in the sense of probability distributions), and (ii) $D_{KL}(p\| q)\neq D_{KL}(q\| p)$, that is the KL-divergence is not symmetric in its arguments.
|
||||
|
||||
|
||||
Variational mean-field theory is a systematic way for constructing such an approximate distribution $q(\bd{s},\bd{\theta})$.
|
||||
The main idea is to choose parameters that minimize the difference between the variational free-energy $F_q(\bd{J},\bd{\theta})$ and the true free-energy $F_p(\bd{J}|\beta)$. We will show in Section \ref{subsubsec:EM} below that the difference between these two free-energies is actually the KL-divergence:
|
||||
\be
|
||||
F_q(\bd{J},\bd{\theta})= F_p(\bd{J},\beta) + D_{KL}(q\|p).
|
||||
\ee
|
||||
This equality, when combined with the non-negativity of the KL-divergence has important consequences. First, it shows the variational free-energy is always larger than the true free-energy, $F_q(\bd{J},\bd{\theta})\geq F_p(\bd{J})$, with equality if and only if $q=p$ (the latter inequality is found in many physics textbooks and is known as the Gibbs inequality). Second, finding the best variational free-energy is equivalent to minimizing the KL divergence $D_{KL}(q\| p)$.
|
||||
|
||||
Armed with these observations, let us now derive a MFT of the Ising model using variational methods. In the simplest MFT of the Ising model, the variational distribution is chosen so that all spins are independent:
|
||||
\begin{equation}
|
||||
q(\bd{s},\bd{\theta}) = \frac{1}{Z_q}\exp\left(\sum_i \theta_i s_i\right) = \prod_i \frac{\mathrm e^{\theta_i s_i}}{2\cosh \theta_i}.
|
||||
\end{equation}
|
||||
In other words, we have chosen a distribution $q$ which factorizes on every lattice site. An important property of this functional form is that we can analytically find a closed-form expression for the variational partition function $Z_q$. This simplicity also comes at a cost: ignoring correlations between spins. These correlations become less and less important in higher dimensions and the MFT ansatz becomes more accurate.
|
||||
|
||||
|
||||
To evaluate the variational free-energy, we make use of Eq.~\eqref{eq:main_var_eq}. First, we need the entropy $H_q$ of the distribution $q$. Since $q$ factorizes over the lattice sites, the entropy separates into a sum of one-body terms
|
||||
\begin{eqnarray}
|
||||
H_q(\bd{\theta}) &=& -\sum_{\{s_i=\pm 1\}}q(\bd{s},\bd{\theta})\log q(\bd{s},\bd{\theta})\nonumber\\
|
||||
&=& -\sum_i q_i\log q_i +(1-q_i)\log(1-q_i),
|
||||
\end{eqnarray}
|
||||
where $q_i = \frac{\mathrm e^{\theta_i}}{2\cosh \theta_i}$ is the probability that spin $s_i$ is in the $+1$ state. Next, we need to evaluate the average of the Ising energy $E(\bd{s},\bd{J})$ with respect to the variational distribution $q$. Although the energy contains bilinear terms, we can still evaluate this average easily, because the spins are independent (uncorrelated) in the $q$ distribution. The mean value of spin $s_i$ in the $q$ distribution, or the on-site magnetization, is given by
|
||||
\begin{equation}
|
||||
\label{eq:MF_self_consistency}
|
||||
m_i = \langle s_i\rangle_q = \sum_{s_i=\pm 1 }s_i \frac{\mathrm e^{\theta_i s_i}}{2\cosh \theta_i} = \tanh(\theta_i).
|
||||
\end{equation}
|
||||
Since the spins are independent, we have
|
||||
\begin{equation}
|
||||
\langle E(\bd{s},\bd{J})\rangle_q = -\frac{1}{2}\sum_{i,j}J_{ij}m_i m_j - \sum_i h_i m_i.
|
||||
\end{equation}
|
||||
The total variational free-energy is
|
||||
\begin{equation*}
|
||||
\beta F_q(\bd{J},\bd{\theta}) = \beta \langle E(\bd{s},\bd{J})\rangle_q - H_q,
|
||||
\end{equation*}
|
||||
and minimizing with respect to the variational parameters ${\bd \theta}$, we obtain
|
||||
\begin{equation}
|
||||
\frac{\partial}{\partial \theta_i} \beta F_q({\bd{J},\bd{\theta}}) = 2\frac{d q_i}{d\theta_i}\left(-\beta\left[\sum_j J_{ij} m_j + h_i\right] + \theta_i \right).
|
||||
\end{equation}
|
||||
Setting this equation to zero, we arrive at
|
||||
\begin{equation}
|
||||
\label{eq:MF_Ising}
|
||||
\theta_i = \beta\sum_j J_{ij}m_j(\theta_j) + h_i.
|
||||
\end{equation}
|
||||
For the special case of a uniform field $h_i = h$ and uniform nearest neighbor couplings $J_{ij} = J$, by symmetry the variational parameters for all the spins are identical, with $\theta_i = \theta$ for all $i$. Then, the mean-field equations reduce to their familiar textbook form, $m = \tanh(\theta)$ and $\theta = \beta(z J m(\theta) + h)$, where $z$ is the coordination number of the lattice (i.e. the number of nearest neighbors).
|
||||
|
||||
Equations~\eqref{eq:MF_self_consistency} and~\eqref{eq:MF_Ising} form a closed system, known as the mean-field equations for the Ising model. To find a solution to these equations, one method is to iterate through and update each $\theta_i$, once at a time, in an asynchronous fashion. To see the relationship of this approach to solving the MFT equations to Expectation Maximization (EM), it is helpful to explicitly spell out the iterative procedure to find the solutions to Eq. ~\eqref{eq:MF_Ising}. We start by initializing our variational parameters to some $\bd{\theta}^{(0)}$ and repeat the following until convergence:
|
||||
\begin{enumerate}
|
||||
\item \emph{Expectation}: Given a set of assignments at iteration $t$, $\bd{\theta}^{(t)}$, calculate the corresponding magnetizations $\bd{m}^{(t)}$ using Eq.~\eqref{eq:MF_self_consistency}
|
||||
\item \emph{Maximization}: Given a set of magnetizations $m_t$, find new assignments $\theta^{(t+1)}$ which minimize the variational free energy $F_q$. From, Eq.~\eqref{eq:MF_Ising}
|
||||
this is just
|
||||
\be
|
||||
\theta_i^{(t+1)}= \beta\sum_j J_{ij}m_j^{(t)} + h_i.
|
||||
\ee
|
||||
\end{enumerate}
|
||||
From these equations, it is clear that we can think of the MFT of the Ising model as an EM-like procedure similar to the one we used for k-means clustering and Gaussian Mixture Models in Sec.~\ref{sec:clustering}.
|
||||
|
||||
|
||||
As if well known in statistical physics, even though MFT is not exact, it can often yield qualitatively and even quantitatively precise predictions (especially in high dimensions). The discrepancy between the true physics and MFT predictions stems from the fact that the variational distribution $q$ we chose does not model the correlations between spins. For instance, it predicts the wrong value for the critical temperature for the $2D$ Ising model. It even erroneously predicts the existence of a phase transition in one dimension at a non-zero temperature. We refer the interested reader to standard textbooks on statistical physics for a detailed analysis of applicability of MFT to the Ising model. However, we emphasize that the failure of any particular variational ansatz does not compromise the power of the approach. In some cases, one can consider changing the variational ansatz to improve the predictive properties of the corresponding variational MFT \cite{yedidia2001idiosyncratic, yedidia2003understanding}. The take-home message is that variational MFT is a powerful tool but one that must be applied and interpreted with care.
|
||||
|
||||
|
||||
|
||||
\subsection{\label{subsubsec:EM}Expectation Maximization (EM)}
|
||||
Ideas along the lines of variational MFT have been independently developed in statistics and imported into machine learning to perform maximum likelihood (ML) estimates. In this section, we explicitly derive the Expectation Maximization (EM) algorithm and demonstrate further its close relation to variational MFT \cite{neal1998view}. We will focus on latent variable models where some of the variables are hidden and cannot be directly observed. This often makes maximum likelihood estimation difficult to implement. EM gets around this difficulty by using an iterative two-step procedure, closely related to variational free-energy based approximation schemes in statistical physics.
|
||||
|
||||
To set the stage for the following discussion, let $\bd{x}$ be the set of visible variables we can directly observe and $\bd{z}$ be the set of latent or hidden variables that we cannot directly observe. Denote the underlying probability distribution from which $\bd{x}$ and $\bd{z}$ are drawn by $p(\bd{z},\bd{x}|\bd{\theta})$, with $\bd{\theta}$ representing all relevant parameters. Given a dataset $\bd{x}$, we wish to find the maximum likelihood estimate of the parameters $\bd{\theta}$ that maximizes the probability of the observed data.
|
||||
|
||||
As in variational MFT, we view $\bd{\theta}$ as variational parameters chosen to maximize the log-likelihood $L(\bd{\theta})=\log p(\bd{x}|\bd{\theta})$. Algorithmically, this can be done by iterating the variational parameters $\bd{\theta}^{(t)}$ in a series of steps ($t=1,2, \dots$) starting from some arbitrary initial value $\bd{\theta}^{(0)}$:
|
||||
\begin{enumerate}
|
||||
\item {\bf Expectation step (E step): }{ Given the known values of observed variable $\bd{x}$ and the current estimate of parameter $\bd\theta_{t-1}$, find the probability distribution of the latent variable $\bd{z}$:}
|
||||
\begin{equation}\label{eq:Estep}
|
||||
q_{t-1}(\bd{z})=p(\bd{z}|\bd{\theta}^{(t-1)},\bd{x})
|
||||
\end{equation}
|
||||
\item {\bf Maximization step (M step):} { Re-estimate the parameter $\bd{\theta}^{(t)}$ to be those with maximum likelihood, assuming $q_{t-1}(\bd z)$ found in the previous step is the true distribution of hidden variable $\bd z$:}
|
||||
\begin{equation}\label{eq:Mstep}
|
||||
\bd{\theta}_{t}=\argmax_{\bd{\theta}}\langle \log p(\bd{z},\bd{x}|\bd{\theta})\rangle_{q_{t-1}}
|
||||
\end{equation}
|
||||
\end{enumerate}
|
||||
It was shown \cite{dempster1977maximum} that {each EM iteration increases the true log-likelihood $L(\bd{\theta})$}, or at worst leaves it unchanged. { In most models, this iteration procedure converges to a \emph{local maximum} of $L(\bd \theta)$.
|
||||
|
||||
|
||||
\begin{figure}[h!]
|
||||
\includegraphics[width=1.0\columnwidth]{variational/EMalg.pdf}
|
||||
\caption{Convergence of EM algorithm. Starting from $\bd{\theta}_t$, E-step (blue) establishes $F_{q}(\bd{\theta}^{(t)})$ which is always a lower bound of $F_p:=\log p(\bd x|\bd \theta)$ (green). M-step (red) is then applied to update the parameter, yielding $\bd{\theta}^{(t+1)}$. The updated parameter $\bd{\theta}^{(t+1)}$ is then used to construct $F_{q}(\bd{\theta}^{(t+1)})$ in the subsequent E-step. M-step is performed again to update the parameter, etc. }
|
||||
\label{fig:EM}
|
||||
\end{figure}
|
||||
|
||||
|
||||
|
||||
To see how EM is actually performed and related to variational MFT, we make use of KL-divergence between two distributions introduced in the last section. Recall that our goal is to maximize the log-likelihood $L(\bd{\theta})=\log p(\bd{x}|\bd{\theta})$. With data $\bd{z}$ missing, we surely cannot just maximize $L(\bd{\theta})$ directly since parameter $\bd\theta$ might couple both $\bd z$ and $\bd x$. EM circumvents this by optimizing another objective function, $F_{q}(\bd{\theta})$, constructed based on estimates of the hidden variable distribution $q(\bd{z}|\bd{x})$. Indeed, the function optimized is none other than the {\it variational free energy} we encountered in the previous section:
|
||||
\begin{eqnarray}
|
||||
F_{q}(\bd\theta) := - \langle \log p(\bd{z},\bd{x}|\bd{\theta})\rangle_{q} -H_q \label{eq:Fptilde},
|
||||
\end{eqnarray}
|
||||
where $H_q$ is the Shannon entropy (defined in Eq.~\eqref{eq:entropy}) of $q(\bd z)$. Once can define the true free-energy $F_p(\bd\theta)$ as the negative log-likelihood of the observed data:
|
||||
\be
|
||||
-F_p(\bd\theta) = L(\bd{\theta}) = \log p(\bd{x}|\bd{\theta}).
|
||||
\ee
|
||||
In the language of statistical physics, $F_p(\bd \theta)$ is the {\it true} free-energy while $F_{q}(\bd\theta)$ is the variational free-energy we would like to minimize (see Table \ref{tbl:vMFT_EM}). Note that we have chosen to employ a physics sign convention here of defining the free-energy as minus log of the partition function. In the ML literature, this minus sign is often omitted \cite{neal1998view} and this can lead to some confusion.
|
||||
Our goal is to choose $\bd{\theta}$ so that our variational free-energy $F_{q}(\bd{\theta})$ is as close to the true free-energy $F_p(\bd{\theta})$ as possible. The difference between these free-energies can be written as
|
||||
\begin{eqnarray}
|
||||
&& F_{q}(\bd\theta) - F_{p}(\bd\theta) \nonumber\\
|
||||
&=& \log p(\bd{x}|\bd{\theta}) -\sum_{\bd z}q(\bd{z}|\bd{x}) \log p(\bd{z},\bd{x}|\bd{\theta}) \nonumber \\
|
||||
+ &&\sum_{\bd z}q(\bd{z}|\bd{x}) \log q(\bd{z}|\bd{x}) \nonumber\\
|
||||
&=&\sum_{\bd z}q(\bd{z}|\bd{x}) \log p(\bd{x}|\bd{\theta}) -\sum_{\bd z}q(\bd{z}|\bd{x}) \log p(\bd{z},\bd{x}|\bd{\theta}) \nonumber\\
|
||||
&& + \sum_{\bd z}q(\bd{z}|\bd{x}) \log q(\bd{z}|\bd{x}) \nonumber\\
|
||||
&=&\! -\sum_{\bd z}q(\bd{z}|\bd{x}) \log \frac{p(\bd{z},\bd{x}|\bd{\theta})}{p(\bd{x}|\bd{\theta})}
|
||||
\!+\!\sum_{\bd z}q(\bd{z}|\bd{x}) \log\tilde p(\bd{z}) \nonumber\\
|
||||
&=& \sum_{\bd z}q(\bd{z}|\bd{x}) \log \frac{q(\bd{z}|\bd{x})}{ p(\bd{z}|\bd{x},\bd{\theta})} \nonumber \\
|
||||
&=& D_{KL}(q(\bd{z}|\bd{x}) \| p(\bd{z}|\bd{x},\bd{\theta}))\ge 0 \nonumber
|
||||
\end{eqnarray}
|
||||
where we have used Bayes' theorem $ p(\bd z| \bd{x},\bd{\theta}) = p (\bd z, \bd x|\bd \theta)/ p (\bd x|\bd \theta)$. Since the KL-divergence is always positive, this shows that the variational free-energy {$F_{q}$ is always an upper bound of the true free-energy $F_p$}. In physics, this result is known as Gibbs' inequality.
|
||||
|
||||
From Eq.~\eqref{eq:Fptilde} and the fact that the the entropy term in Eq.~\eqref{eq:Fptilde} does not depend on $\bd\theta$, we can immediately see that the maximization step (M-step) in Eq.~\eqref{eq:Mstep} is equivalent to minimizing the variational free-energy $F_{q}(\bd\theta)$. Surprisingly, the expectation step (E-step) can also viewed as the optimization of this variational free-energy. Concretely, one can show that the distribution of hidden variables $\bd{z}$ given the observed variable $\bd x$ and the current estimate of parameter $\bd\theta$, Eq.~(\ref{eq:Estep}), is the \emph{unique} probability $q(\bd z)$ that minimizes $F_{q}(\bd\theta)$ (now seen as a functional of $q$). This can be proved by taking the functional derivative of Eq.~\eqref{eq:Fptilde}, plus a Lagrange multiplier that encodes $\sum_{\bd{z}} q(\bd z)=1$, with respect to $q(\bd z)$. Summing things up, we can re-write EM in the following form \cite{neal1998view}:
|
||||
\begin{enumerate}
|
||||
\item \emph{Expectation step:} Construct the approximating probability distribution {of unobserved $\bd z$ given the values of observed variable $\bd x$ and parameter estimate $\bd\theta^{(t-1)}$}:
|
||||
\begin{equation}\label{eq:EM_Estep}
|
||||
q_{t-1}(\bd{z})=\argmin_{q}F_{q}(\bd{\theta}^{(t-1)})
|
||||
\end{equation}
|
||||
\item \emph{Maximization step:} {Fix $q$, }update the variational parameters:
|
||||
\begin{equation}\label{eq:EM_Mstep}
|
||||
\bd{\theta}^{(t)}=\argmax_{\bd{\theta}} -F_{q_{t-1}}(\bd{\theta}).
|
||||
\end{equation}
|
||||
\end{enumerate}
|
||||
|
||||
|
||||
To recapitulate, EM implements ML estimation even with missing or hidden variables through optimizing a lower bound of the true log-likelihood. In statistical physics, this is reminiscent of optimizing a variational free-energy which is a lower bound of true free-energy due to Gibbs inequality. In Fig.~\ref{fig:EM}, we show pictorially how EM works. The E-step can be seen as representing the unobserved variable $\bd z$ by a probability distribution $q(\bd z)$. This probability is used to construct an alternative objective function $-F_{q}(\bd\theta)$, which is then maximized with respect to $\bd\theta$ in the M-step. By construction, maximizing the negative variational free-energy is equivalent to doing ML estimation on the joint data (i.e.~both observed and unobserved). The name ``M-step'' is intuitive since the parameters $\bd\theta$ are found by maximizing $-F_{q}(\bd\theta)$. The name ``E-step'' comes from the fact that one usually doesn't need to construct the probability of missing datas explicitly, but rather need only compute the ``expected" sufficient statistics over these data, cf.~Fig.~\ref{fig:EM}.
|
||||
|
||||
|
||||
|
||||
On the practical side, EM has been demonstrated to be extremely useful in parameter estimation, particularly in hidden Markov models and Bayesian networks (see, for example, \cite{wainwright2008graphical,barber2012bayesian}). One of the striking advantages of EM is that it is conceptually simple and easy to implement (see \href{https://physics.bu.edu/~pankajm/MLnotebooks.html}{Notebook 16}). In many cases, implementation of EM is guaranteed to increase the likelihood monotonically, which could be a perk during debugging. For readers interested in an overview on applications of EM, we recommend \cite{do2008expectation}.
|
||||
|
||||
Finally for advanced readers familiar with the physics of disordered systems, we note that it is possible to construct a one-to-one dictionary between EM for latent variable models and the MFT of spin systems with quenched disorder. In a disordered spin systems, the Ising couplings $\bd{J}$ are commonly taken to be quenched random variables drawn from some underlying probability distribution. In the EM procedure, the quenched disorder is provided by the observed data points $\bd{x}$ which are drawn from some underlying probability distribution that characterizes the data. The spins $\bd{s}$ are like the hidden or latent variables $\bd{z}$. Similar analogues can be found for all the variational MFT quantities (see Table \ref{tbl:vMFT_EM}). This striking correspondence offers a glimpse into the deep connection between statistical mechanics and unsupervised latent variable models -- a connection that we will repeatedly exploit to gain more intuition for the energy-based unsupervised models considered in the next few chapters.
|
||||
|
||||
\begin{table}[t!]
|
||||
\begin{tabular}{*{3}{|p{0.5\columnwidth}}}
|
||||
\hline
|
||||
\emph{\bf{statistical physics}} & \emph{\bf{Variational EM}} \\
|
||||
\hline
|
||||
spins/d.o.f.: $\bd{s}$ & hidden/latent variables $\bd{z}$ \\
|
||||
\hline
|
||||
couplings /quenched disorder: $\bd{J}$ & data observations: $\bd{x}$ \\
|
||||
\hline
|
||||
Boltzmann factor $e^{-\beta E(\bd{s},\bd{J})}$ & Complete probability: $p({\bd x},{\bd z}|\bd{\theta})$ \\
|
||||
\hline
|
||||
partition function: $Z(\bd{J})$ & marginal likelihood $p(\bd{x}|\theta)$ \\
|
||||
\hline
|
||||
energy: $\beta E(\bd{s},\bd{J})$ & negative log-complete data likelihood: $-\log p({\bd x},{\bd z}|\bd{\theta},m)$ \\
|
||||
\hline
|
||||
free energy: $\beta F_p(\bd{J}|\beta)$ & negative log-marginal likelihood: $-\log p(\bd{x}|m)$ \\
|
||||
\hline
|
||||
|
||||
variational distribution: $q(\bd{s})$ & variational distribution: $q(\bd{z}|\bd{x})$ \\
|
||||
\hline
|
||||
variational free-energy: $F_q(\bd{J}, \bd{\theta}) $ & variational free-energy:
|
||||
$F_{q}(\bd{x}, \bd\theta)$ \\
|
||||
\hline
|
||||
\end{tabular}
|
||||
\caption{Analogy between quantities in statistical physics and variational EM.}
|
||||
\label{tbl:vMFT_EM}
|
||||
\end{table}
|
||||
|
||||
|
||||
%\subsection{Graphical Models, Bp, and variational methods}
|
||||
|
||||
|
||||
|
||||
@@ -134,3 +377,462 @@ process then follows the same structure as the VMC method.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
\section{Energy Based Models: Maximum Entropy (MaxEnt) Principle, Generative models, and Boltzmann Learning}
|
||||
\label{sec:EnergyI}
|
||||
|
||||
Most of the models discussed in the previous sections (e.g. linear and logistic regression, ensemble models, and supervised neural networks) are \emph{discriminative} -- they are designed to perceive differences between groups or categories of data. For example, recognizing differences between images of cats and images of dogs allows a discriminative model to label an image as ``cat'' or ``dog''. Discriminative models form the core techniques of most supervised learning methods. However, discriminative methods have several limitations. First, like all supervised learning methods, they require labeled data. Second, there are tasks that discriminative approaches simply cannot accomplish, such as drawing new examples from an unknown probability distribution. A model that can learn to represent and sample from a probability distribution is called \emph{generative}. For example, a generative model for images would learn to draw new examples of cats and dogs given a dataset of images of cat and dog. Similarly, given samples generated from one phase of an Ising model we may want to generate new samples from that phase. Such tasks are clearly beyond the scope of discriminative models like the ensemble models and DNNs discussed so far in the review. Instead, we must turn to a new class of machine learning methods.
|
||||
|
||||
The goal of this section is to introduce the reader to \emph{energy-based} generative models. As we will see, energy-based models are closely related to the kinds of models commonly encountered in statistical physics. For this reason, we will draw upon many techniques that have their origin in statistical mechanics (e.g. Monte-Carlo methods). The section starts with a brief overview of generative models, highlighting the similarities and differences with the supervised learning methods encountered in earlier sections. Next, we introduce perhaps the simplest kind of generative model -- Maximum Entropy (MaxEnt) models. MaxEnt models have no latent (or hidden) variables, making them ideal for introducing the key concepts and tools that underlie energy-based generative models. We then present an extended discussion of how to train energy-based models. Much of this discussion will be also be applicable to more complicated energy-based models such as Restricted Boltzmann Machines (RBMs) and the deep models discussed in the next section.
|
||||
\subsection{An overview of energy-based generative models}
|
||||
|
||||
Generative models are machine learning methods that learn to generate new examples similar to those found in a training dataset. The core idea of most generative models is to learn a parametric model for the probability distribution from which the data was drawn. Once we have learned a model, we can generate new examples by sampling from the learned generative model (see Figure~\ref{fig:MNIST-generative}). As in statistical physics, this sampling is often done using Markov Chain Monte Carlo (MCMC) methods. A review of MCMC methods is beyond the scope of this discussion: for a concise and beautiful introduction to MCMC-inspired methods that bridges both statistical physics and ML the reader is encouraged to consult Chapters 29-32 of David MacKay's book \cite{mackay2003information} as well as the review by Michael I. Jordan and others~\cite{andrieu2003introduction}.
|
||||
|
||||
|
||||
The added complexity of learning models directly from samples introduces many of the same fundamental tensions we encountered when discussing discriminative models. The ability to generate new examples requires models to be able to ``generalize'' beyond the examples they have been trained on, that is to generate new samples that are not samples of the training set. The models must be expressive enough to capture the complex correlations present in the underlying data distribution, but the amount of data we have is finite which can give rise to overfitting.
|
||||
|
||||
|
||||
|
||||
In practice, most generative models that are used in machine learning flexible enough that, with a sufficient number of parameters, they can approximate any probability distribution. For this reason, there are three axes on which we can differentiate classes of generative models:
|
||||
\begin{itemize}
|
||||
\item The first axis is how easy the model is to train -- both in terms of computational time and the complexity of writing code for the algorithm.
|
||||
\item The second axis is how well the model generalizes from the training set to the test set.
|
||||
\item The third axis is which characteristics of the data distribution the model is capable of and focuses on capturing.
|
||||
\end{itemize}
|
||||
All generative models must balance these competing requirements and generative models differ in the tradeoffs they choose. Simpler models capture less structure about the underlying distributions but are often easier to train. More complicated models can capture this structure but may overfit to the training data.
|
||||
|
||||
|
||||
|
||||
One of the fundamental reasons that energy-based models have been less widely-employed than their discriminative counterparts is that the training procedure for these models differs significantly from those for supervised neural networks models. Though both employ gradient-descent based procedures for minimizing a cost function (one common choice for generative models is the negative log-likelihood function), energy-based models do not use backpropagation (see Sec.~\ref{subsec:backprop}) and automatic differentiation for computing gradients. Rather, one must turn to ideas inspired by MCMC based methods in physics and statistics that sometimes go under the name ``Boltzmann Learning'' (discussed below). As a result, training energy-based models requires additional tools that are not immediately available in packages such as PyTorch and TensorFlow.
|
||||
|
||||
|
||||
A new open-source package -- {\it Paysage} -- that is built on top of PyTorch bridges this gap by providing the toolset for training energy-based models (Paysage is maintained by Unlearn.AI -- a company affiliated with two of the authors (CKF and PM)). Paysage makes it easy to quickly code and deploy energy-based models such as Restricted Boltzmann Machines (RBMs) and Stacked RBMs-- a ``deep'' unsupervised model. The package includes unpublished training methods that significantly improve the training performance, can be used for various datatypes, and can be employed on GPUs. We make use of this package extensively in the next two sections and the accompanying Python notebooks. For example, Figure \ref{fig:MNIST-generative} (and the accompanying \href{https://physics.bu.edu/~pankajm/MLnotebooks.html}{notebook}) show how the Paysage package can be used to quickly code and train a variety of energy-based models on the MNIST handwritten digit dataset.
|
||||
|
||||
|
||||
|
||||
|
||||
Finally, we note that generative models at their most basic level are complex parameterizations of the probability distribution the data is drawn from. For this reason, generative models can do much more than just generate new examples. They can be used to perform a multitude of other tasks that require sampling from a complex probability distribution including ``de-noising'', filling in missing data, and even discrimination \cite{hinton2012practical}. The versatility of generative models is one of the major appeals of these unsupervised learning methods.
|
||||
|
||||
\begin{figure}[t]
|
||||
\includegraphics[width=0.9\columnwidth]{Energy/mnist-RBM.png}
|
||||
\caption{Examples of handwritten digits (``reconstructions'') generated using various energy-based models using the powerful \emph{Paysage} package for unsupervised learning. Examples from top to bottom are: the original MNIST database, an RBM with Gaussian units which is equivalent to a Hopfield Model, a Restricted Boltzmann Machine (RBM), a RBM with an $L_1$ penalty for regularization, and a Deep Boltzmann Machine (DBM) with $3$ layers. All models have $200$ hidden units. See Sec. \ref{sec:EnergyII} and corresponding notebook for details}
|
||||
\label{fig:MNIST-generative}
|
||||
\end{figure}
|
||||
|
||||
|
||||
\subsection{Maximum entropy models: the simplest energy-based generative models}
|
||||
|
||||
Maximum Entropy (MaxEnt) models are one of the simplest classes of energy-based generative models. MaxEnt models have their origin in a series of beautiful papers by Jaynes that reformulated statistical mechanics in information theoretic terms \cite{jaynes1957information, jaynes1957informationII}. Recently, the flood of new, large scale datasets has resulted in a resurgence of interest in MaxEnt models in many fields including physics (especially biological physics), computational neuroscience, and ecology \cite{schneidman2006weak, weigt2009identification, elith2011statistical}. MaxEnt models are often presented as the class of generative models that make the least assumptions about the underlying data. However, as we have tried to emphasize throughout the review, all ML and statistical models require assumptions, and MaxEnt models are no different. Overlooking this can sometimes lead to misleading conclusions, and it is important to be cognizant of these implicit assumptions \cite{schwab2014zipf, aitchison2016zipf}.
|
||||
|
||||
|
||||
|
||||
\subsubsection{MaxEnt models in statistical mechanics}
|
||||
|
||||
MaxEnt models were introduced by E.~T.~Jaynes in a two-part paper in 1957 entitled ``Information theory and statistical mechanics'' \cite{jaynes1957information, jaynes1957informationII}. In these incredible papers, Jaynes showed that it was possible to re-derive the Boltzmann distribution (and the idea of generalized ensembles) entirely from information theoretic arguments. Quoting from the abstract, Jaynes considered ``statistical mechanics as a form of statistical inference rather than as a physical theory'' (portending the close connection between statistical physics and machine learning). Jaynes showed that the Boltzmann distribution could be viewed as resulting from a statistical inference procedure for learning probability distributions describing physical systems where one only has partial information about the system (usually the average energy).
|
||||
|
||||
|
||||
The key quantity in MaxEnt models is the information theoretic, or Shannon, entropy, a concept introduced by Shannon in his landmark treatise on information theory \cite{shannon1949communication}. The Shannon entropy quantifies the statistical uncertainty one has about the value of a random variable ${\bf x}$ drawn from a probability distribution $p({\bf x})$. The Shannon entropy of the distribution is defined as
|
||||
\be
|
||||
S_p = - \mathrm{Tr}_{\bf x} p({\bf x}) \log p({\bf x})
|
||||
\ee
|
||||
where the trace is a sum/integral over all possible values a variable can take. Jaynes showed that the Boltzmann distribution follows from the Principle of Maximum Entropy. A physical system should be described by the probability distribution with the largest entropy subject to certain constraints (often provided by measuring the average value of conserved, extensive quantities such as the energy, particle number, etc.) The principle uniquely specifies a procedure for parameterizing the functional form of the probability distribution. Once we have specified and learned this form we can, of course, generate new examples by sampling this distribution.
|
||||
|
||||
|
||||
Let us illustrate how this works in more detail. Suppose that we have chosen a set of functions $\{f_i({\bf x})\}$ whose average value we want to fix to some observed values $\langle f_i \rangle_\mathrm{obs}$. The Principle of Maximum Entropy states that we should choose the distribution $p({\bf x})$ with the largest uncertainty (i.e.~largest Shannon entropy $S_p$), subject to the constraints that the model averages match the observed averages:
|
||||
\be
|
||||
\langle f_i \rangle_\mathrm{model} := \int \mathrm{d} {\bf x} f_i({\bf x}) p({\bf x}) = \langle f_i \rangle_\mathrm{obs}
|
||||
\ee
|
||||
|
||||
We can formulate the Principle of Maximum Entropy as an optimization problem using the method of Lagrange multipliers by minimizing:
|
||||
\begin{align}
|
||||
\mathcal{L}[p] &= -S_p + \sum_i \lambda_i \left(\langle f_i \rangle_\mathrm{obs} - \int d {\bf x} f_i({\bf x}) p({\bf x})\right) \nonumber \\
|
||||
&+ \gamma \left(1 - \int \mathrm{d} {\bf x} p({\bf x}) \right) \nonumber,
|
||||
\end{align}
|
||||
where the first set of constraints enforce the requirement for the averages and the last constraint enforces the normalization that the trace over the probability distribution equals one. We can solve for $p({\bf x})$ by taking the functional derivative and setting it to zero
|
||||
\begin{align}
|
||||
0=\frac{\delta \mathcal{L}}{ \delta p} = (\log p({\bf x})) + 1) -\sum_i \lambda_i f_i({\bf x}) - \gamma \nonumber.
|
||||
\end{align}
|
||||
The general form of the maximum entropy distribution is then given by
|
||||
\begin{equation}
|
||||
p({\bf x}) = \frac{1}{Z}\; \mathrm e^{\;\sum_i \lambda_i f_i({\bf x})}
|
||||
\label{eq:boltzmann-dist}
|
||||
\end{equation}
|
||||
where $Z({\lambda_i}) = \int \mathrm{d} {\bf x}\mathrm e^{\sum_i \lambda_i f_i({\bf x})}$ is the partition function.
|
||||
|
||||
The maximum entropy distribution is clearly just the usual Boltzmann distribution with energy $E({\bf x}) = - \sum_i \lambda_i f_i({\bf x})$. The values of the Lagrange multipliers are chosen to match the observed averages for the set of functions $\{f_i({\bf x})\}$ whose average value is being fixed:
|
||||
\begin{equation}
|
||||
\langle f_i \rangle_\mathrm{model} = \int \mathrm{d} {\bf x} p({\bf x}) f_i({\bf x})= \frac{\partial \log Z}{\partial \lambda_i} = \langle f_i \rangle_\mathrm{obs}.
|
||||
\label{eq:moment-constraint}
|
||||
\end{equation}
|
||||
In other words, the parameters of the distribution can be chosen such that
|
||||
\be
|
||||
\partial_{\lambda_i} \log Z = \langle f_i \rangle_\mathrm{data}.
|
||||
\ee
|
||||
|
||||
To gain more intuition for the MaxEnt distribution, it is helpful to relate the Lagrange multipliers to the familiar thermodynamic quantities we use to describe physical systems \cite{jaynes1957information}. Our $\mathbf{x}$ denotes the microscopic state of the system, i.e. the MaxEnt distribution is a probability distribution over microscopic states. However, in thermodynamics we only have access to average quantities. If we know only the average energy $ \langle E(\mathbf{x}) \rangle_\mathrm{obs}$, the MaxEnt procedure tells us to maximize the entropy subject to the average energy constraint. This yields
|
||||
\be
|
||||
p({\mathbf x})= \frac{1}{Z}\; \mathrm e^{- \beta E(\mathbf x)}.
|
||||
\ee
|
||||
where we have identified the Lagrange multiplier conjugate to the energy $\lambda_1=-\beta= 1/k_BT$ with the (negative) inverse temperature. Now, suppose we also constrain the particle number $\langle N(x) \rangle_\mathrm{obs}$. Then, an almost identical calculation yields a MaxEnt distribution of the functional form
|
||||
\be
|
||||
p({\mathbf x})= \frac{1}{Z}\; \mathrm e^{-\beta( E({\mathbf x}) -\mu N(\mathbf{x}))},
|
||||
\ee
|
||||
where we have rewritten our Lagrange multipliers in the familiar thermodynamic notation $\lambda_1= -\beta$ and $\lambda_2 = \mu/\beta$. Since this is just the Boltzmann distribution, we can also relate the partition function in our MaxEnt model to the thermodynamic free-energy via $F= -\beta^{-1} \log Z$. The choice of which quantities to constrain is equivalent to working in different thermodynamic ensembles.
|
||||
|
||||
|
||||
|
||||
\subsubsection{From statistical mechanics to machine learning}
|
||||
|
||||
|
||||
The MaxEnt idea also provides a general procedure for learning a generative model from data. The key difference between MaxEnt models in (theoretical) physics and ML is that in ML we have no direct access to observed values $\langle f_i \rangle_\mathrm{obs}$. Instead, these averaged must be directly estimated from data (samples). To denote this difference, we will call empirical averages calculated from data as $\langle f_i \rangle_\mathrm{data}$. We can think of MaxEnt as a statistical inference procedure simply by replacing $\langle f_i \rangle_\mathrm{obs}$ by $\langle f_i \rangle_\mathrm{data}$ above.
|
||||
|
||||
This subtle change has important implications for training MaxEnt models. First, since we do not know these averages exactly, but must estimate from the data, our training procedures must be careful not to overfit to the observations (our samples might not be reflective of the true values of these statistics). Second, the averages of certain functions $f_i$ are easier to estimate from limited data than others. This is often an important consideration when formulating which MaxEnt model to fit to the data. Finally, we note that unlike in physics where conservation laws often guide the functions $f_i$ whose averages we hold fix, ML offers no comparable guide for how to choose the $f_i$ we care about. For these reasons, choosing the $\{ f_i \}$ is often far from straightforward.
|
||||
As a final point, we note that here we have presented a physics-based perspective for justifying the MaxEnt procedure. We mention in passing that the MaxEnt in ML is also closely related to ideas from Bayesian inference \cite{jaynes1968prior, jaynes2003probability} and this latter point of view is more common in discussion of MaxEnt in the statistics and ML literature.
|
||||
|
||||
|
||||
|
||||
\subsubsection{Generalized Ising Models from MaxEnt}
|
||||
|
||||
The form of a MaxEnt model is completely specified once we choose the averages $\{f_i\}$ we wish to constrain. One common choice often used in MaxEnt modeling is to constrain the first two moments of a distribution. When our random variables $\bf{x}$ are continuous, the corresponding MaxEnt distribution is a multi-dimensional Gaussian. If the $\bf{x}$ are binary (discrete), then the corresponding MaxEnt distribution is a generalized Ising (Potts) model with all-to-all couplings.
|
||||
|
||||
|
||||
To see this, consider a random variable $\bf{x}$ with first and second moments $\langle x_i \rangle_\mathrm{data}$ and $\langle x_i x_j \rangle_\mathrm{data}$, respectively. According to the Principle of Maximum Entropy, we should choose to model this variable using a Boltzmann distribution with constraints on the first and second moments. Let $a_i$ be the Lagrange multiplier associated with $\langle x_i \rangle_\mathrm{data}$ and $J_{ij} / 2$ be the Lagrange multiplier associated with $\langle x_i x_j \rangle_\mathrm{data}$. Using Eq.~\eqref{eq:moment-constraint}, it is easy to verify that the energy function
|
||||
\begin{equation}
|
||||
E({\bf x}) = -\sum_i a_i x_i - \frac{1}{2} \sum_{ij} J_{ij} x_i x_j
|
||||
\label{eq:pairwise-maxent}
|
||||
\end{equation}
|
||||
satisfies the above constraints.
|
||||
|
||||
|
||||
|
||||
Partition functions for maximum entropy models are often intractable to compute. Therefore, it is helpful to consider two special cases where ${\bf x}$ has different support (different kinds of data). First, consider the case that the random variables ${\bf x} \in \mathbb{R}^n$ are real numbers. In this case we can compute the partition function directly:
|
||||
\begin{align}
|
||||
Z
|
||||
&= \int \mathrm{d} {\bf x} \mathrm e^{{\bf a}^T {\bf x} + \frac{1}{2} {\bf x}^T J {\bf x} }
|
||||
= \sqrt{(2 \pi)^n \text{det} J^{-1}} \mathrm e^{-\frac{1}{2} {\bf a}^T J^{-1} {\bf a}}.
|
||||
\end{align}
|
||||
The resulting probability density function is,
|
||||
\begin{align}
|
||||
p({\bf x}) &= Z^{-1} \mathrm e^{ -E({\bf x})} \nonumber \\
|
||||
&= \frac{1}{ \sqrt{(2 \pi)^n \text{det} J^{-1}}} \mathrm e^{\frac{1}{2} {\bf a}^T J^{-1} {\bf a} + {\bf a}^T {\bf x} + \frac{1}{2} {\bf x}^T J {\bf x} } \nonumber \\
|
||||
&= \frac{1}{ \sqrt{(2 \pi)^n \text{det} \Sigma }} \mathrm e^{-\frac{1}{2} ({\bf x} - {\bf \mu})^T \Sigma^{-1} ( {\bf x} - {\bf \mu}) },
|
||||
\label{eq:gaussian-pdf}
|
||||
\end{align}
|
||||
where ${\bf \mu} = -J^{-1} {\bf a}$ and $\Sigma = -J^{-1}$. This, of course, is the normalized, multi-dimensional Gaussian distribution.
|
||||
|
||||
Second, consider the case that the random variable ${\bf x}$ is binary with $x_i \in \{-1, +1\}$. The energy function takes the same form as Eq.~\eqref{eq:pairwise-maxent}, but the partition function can no longer be computed in a closed form. This model is known as the Ising model in the physics literature, and is often called a Markov Random Field in the machine learning literature. It is well known to physicists that calculating the partition function for the Ising Model is intractable. For this reason, the best we can do is estimate it using numerical techniques such MCMC methods or approximate methods like variational MFT methods, see Sec.~\ref{sec:varl_MFT}. Finally, we note that in ML it is common to use binary variables which take on values in $x_i \in \{0,1\}$ rather than $\{\pm 1\}$. This can sometimes be a source of confusion when translating between ML and physics literatures and can lead to confusion when using ML packages for physics problems.
|
||||
|
||||
\subsection{Cost functions for training energy-based models}
|
||||
|
||||
The MaxEnt procedure gives us a way of parameterizing an energy-based generative model. For any energy-based generative model, the energy function $E({\bf x}, \{\theta_i\})$ depends on some parameters $\theta_i$ -- couplings in the language of statistical physics -- that must be inferred directly from the data. For example, for the MaxEnt models the $\{ \theta_i\}$ are just the Lagrange multipliers $\{ \lambda_i\}$ introduced in the last section. The goal of the training procedures is to use the available training data to fit these parameters.
|
||||
|
||||
Like in many other ML techniques, the way we will fit these couplings is by minimizing a cost function using stochastic gradient descent (cf.~Sec.~\ref{sec:gradient_descent}). Such a procedure naturally separates into two parts: choosing an appropriate cost function and calculating the gradient of the cost function with respect to the model parameters. Formulating a cost function for generative models is a little bit trickier than for supervised, discriminative models. The objective of discriminative models is straightforward -- predict the label from the features. However, what we mean by a ``good'' generative model is much harder to define using a cost function. We would like the model to generate examples similar to those we find in the training dataset. However, we would also like the model to be able to generalize -- we do not want the model to reproduce ``spurious details'' that are particular to the training dataset. Unlike for discriminative models, there is no straightforward idea like cross-validation on the data labels that neatly addresses this issue. For this reason, formulating cost functions for generative models is subtle and an important and interesting open area of research.
|
||||
|
||||
Calculating the gradients of energy-based models also turns out to be different than for discriminative models such as deep neural networks. Rather than rely on automatic differentiation techniques and backpropagation (see Sec.~\ref{subsec:backprop}), calculating the gradient requires drawing on intuitions from MCMC based methods. Below, we provide an in-depth discussion of Boltzmann learning for energy-based generative models, focusing on MaxEnt models. We put the emphasis on training procedures that generalize to more complicated generative models with latent variables such as RBMs discussed in the next section. Therefore, we largely ignore the incredibly rich physics-based literature on fitting Ising-like MaxEnt models (see the recent reviews \cite{nguyen2017inverse,baldassi2018inverse} and references therein).
|
||||
|
||||
|
||||
|
||||
|
||||
\subsubsection{Maximum likelihood }
|
||||
|
||||
|
||||
By far the most common approach used for training a generative model is to maximize the log-likelihood of the training data set. Recall, that the log-likelihood characterizes the log-probability of generating the observed data using our generative model. By choosing the negative log-likelihood as the cost function, the learning procedure tries to find parameters that maximize the probability of the data. This cost function is intuitive and, for this reason, has been the work-horse of most generative modeling. However, we note that the Maximum Likelihood estimation (MLE) procedure has some important limitations that we will return to in Section \ref{sec:vae}.
|
||||
|
||||
In what follows, we employ a general notation that is applicable to all energy-based models, not just the MaxEnt models introduced above. The reason for this is that much of this discussion does not rely on the specific form of the energy function but only on the fact that our generative model takes a Boltzmann form. We denote the generative model by the probability distribution $p_\theta( {\bf x})$ and its corresponding partition function by $\log{Z(\{\theta_i\})}$. In MLE, the parameters of the model are fit by maximizing the log-likelihood:
|
||||
\bea
|
||||
\mathcal{L}(\{ \theta_i \}) &:=& \langle \log \left( p_\theta( {\bf x}) \right) \rangle_\mathrm{data} \nonumber \\
|
||||
&=& -\langle E({\bf x}; \{ \theta_i\}) \rangle_\mathrm{data} - \log{Z(\{\theta_i\})},
|
||||
\label{Eq:E1_likelihood}
|
||||
\eea
|
||||
where we have set $\beta=1$. In writing this expression we made use of the fact that our generative distribution is of the Boltzmann form and the fact that our partition function does not depend on the data:
|
||||
\be
|
||||
\langle \log Z(\{\theta_i\}) \rangle_\mathrm{data}=\log{Z(\{\theta_i\})}.
|
||||
\ee
|
||||
|
||||
|
||||
|
||||
\subsubsection{Regularization}
|
||||
|
||||
Just as for discriminative models like a linear and logistic regression, it is common to supplement the log-likelihood with additional regularization terms that prevent overfitting (see Secs.~\ref{sec:lin_reg} and \ref{sec:log_reg}). Instead of minimizing the negative log-likelihood, one minimizes a cost function of the form
|
||||
\be
|
||||
-\mathcal{L}(\{ \theta_i \}) + {E}_\mathrm{reg}(\{\theta_i\}),
|
||||
\ee
|
||||
where ${E}_\mathrm{reg}(\{\theta_i\})$ is an additional regularization term that prevents overfitting. From a Bayesian perspective, this new term can be viewed as encoding a (negative) log-prior on model parameters and performing a maximum-a-posteriori (MAP) estimate instead of a MLE (see corresponding discussion in Sec.~\ref{sec:lin_reg}).
|
||||
|
||||
|
||||
|
||||
As we saw with regression, different forms of regularization give rise to different kinds of properties. A common choice for the regularization function are the sums of the $L_1$ or $L_2$ norms of the parameters
|
||||
\be
|
||||
E_\mathrm{reg}(\{\theta_i\})= \Lambda \sum_i |\theta_i|^\alpha, \, \alpha=1,2
|
||||
\ee
|
||||
with $\Lambda$ denoting parameters that control the regularization strength. For $\Lambda=0$, there is no regularization and we are simply performing MLE. In contrast, a choice of large $\Lambda$ will force many parameters to be close to or exactly zero. Just as in regression, an $L1$ penalty enforces sparsity, with many of the $\theta_i$ are set to zero, and $L2$ regularization shrinks the size of the parameters towards zero.
|
||||
|
||||
|
||||
One challenge of generative models is that it is often difficult to choose the regularization strength $\Lambda$. Recall, that for linear and logistic regression $\Lambda$ is chosen to maximize the out-of-sample performance on a validation dataset (i.e.~cross-validation). However, for generative models our data are usually unlabeled. For this reason, choosing a regularization strength is more subtle and there exists no universal procedure for choosing $\Lambda$. One common strategy is to divide the data into a training set and a validation set and monitor a summary statistic such as the log-likelihood, energy distance \cite{szekely2003statistics}, or variational free-energy of the generative model on the training and validation sets (the variational free-energy was discussed extensively in Sec.~\ref{sec:varl_MFT} ) \cite{hinton2012practical}. If the gap between the training and validation datasets starts growing, one is probably overfitting the model even if the log-likelihood of the training dataset is still increasing. This also gives a procedure for ``early stopping'' -- a regularization procedure we introduced in the context of discriminative models. In practice, when using such regularizers it is important to try many different values of $\Lambda$ and then try to use a proxy statistic for overfitting to evaluate the optimal choice of $\Lambda$.
|
||||
|
||||
|
||||
|
||||
|
||||
\subsection{Computing gradients }
|
||||
\label{subsec:energy_gradients}
|
||||
|
||||
We still need to specify a procedure for minimizing the cost function. One powerful and common choice that is widely employed when training energy-based models is stochastic gradient descent (SGD) (see Sec.~\ref{sec:gradient_descent}). Performing MLE using SGD requires calculating the gradient of the log-likelihood Eq.~\eqref{Eq:E1_likelihood} with respect to the parameters $\theta_i$. To simplify notation and gain intuition, it is helpful to define ``operators'' $O_i(\mathbf{x})$, conjugate to the parameters $\theta_i$
|
||||
\be
|
||||
O_i(\mathbf{x}) := {\partial E({\bf x}; \theta_i) \over \partial \theta_i}.
|
||||
\ee
|
||||
Since the partition function is just the cumulant generative function for the Boltzmann distribution, we know that the usual statistical mechanics relationships between
|
||||
expectation values and derivatives of the log-partition function hold:
|
||||
\be
|
||||
\langle O_i (\mathbf{x}) \rangle_\mathrm{model}= \mathrm{Tr}_{\mathbf x} p_\theta(\mathbf{x}) O_i (\mathbf{x}) = -{ \partial \log{Z(\{\theta_i\})} \over \partial \theta_i}.
|
||||
\ee
|
||||
In terms of the operators $\{ O_i(\mathbf{x})\}$, the gradient of Eq.~\eqref{Eq:E1_likelihood} takes the form \cite{ackley1987learning}
|
||||
\bea
|
||||
-{\partial \mathcal{L}(\{\theta_i \}) \over \partial \theta_i}&=& \Big\langle { \partial E({\bf x}; \theta_i) \over \partial \theta_i} \Big\rangle_\mathrm{data}+ {\partial \log{Z(\{\theta_i\})} \over \partial \theta_i} \nonumber \\
|
||||
&=& \langle O_i (\mathbf{x}) \rangle_\mathrm{data}-\langle O_i (\mathbf{x}) \rangle_\mathrm{model}.
|
||||
\label{eq:Pos-neg-energy}
|
||||
\eea
|
||||
|
||||
|
||||
|
||||
These equations have a simple and beautiful interpretation. The gradient of the log-likelihood with respect to a model parameter is a difference of moments -- one calculated directly from the data and one calculated from our model using the current model parameters. The data-dependent term is known as the \emph{positive phase} of the gradient and the model-dependent term is known as the \emph{negative phase} of the gradient. This derivation also gives an intuitive explanation for likelihood-based training procedures. The gradient acts on the model to lower the energy of configurations that are near observed data points while raising the energy of configurations that are far from observed data points. Finally, we note that all information about the data only enters the training procedure through the expectations $\langle O_i (\mathbf{x}) \rangle_\mathrm{data}$ and our generative model is blind to information beyond what is contained in these expectations.
|
||||
|
||||
|
||||
To use SGD, we must still calculate the expectation values that appear in Eq.~\eqref{eq:Pos-neg-energy}. The positive phase of the gradient -- the expectation values with respect to the data -- can be easily calculated using samples from the training dataset. However, the negative phase -- the expectation values with respect to the model -- are generally much more difficult to compute. We will see that in almost all cases, we will have to resort to either numerical or approximate methods. The fundamental reason for this is that it is impossible to calculate the partition function exactly for most interesting models in both physics and ML.
|
||||
|
||||
|
||||
There are exceptional cases in which we can calculate expectation values analytically. When this happens, the generative model is said to have a \emph{Tractable Likelihood}. One example of a generative model with a Tractable Likelihood is the Gaussian MaxEnt model for real valued data discussed in Eq.~\eqref{eq:gaussian-pdf}. The parameters/Lagrange multipliers for this model are the local fields $\mathbf{a}$ and the pairwise coupling matrix $J$. In this case, the usual manipulations involving Gaussian integrals allow us to exactly find the parameters $\mu=-J^{-1}\mathbf{a}$ and $\Sigma=-J^{-1}$, yielding the familiar expressions ${\bf \mu} = \langle {\bf x} \rangle_\mathrm{data}$ and $\Sigma = \langle ({\bf x} - \langle {\bf x} \rangle_\mathrm{data}) ({\bf x} - \langle {\bf x} \rangle_\mathrm{data})^T \rangle_\mathrm{data}$. These are the standard estimates for the sample mean and covariance matrix. Converting back to the Lagrange multipliers yields
|
||||
\be
|
||||
J = - \langle ({\bf x} - \langle {\bf x} \rangle_\mathrm{data}) ({\bf x} - \langle {\bf x} \rangle_\mathrm{data})^T \rangle_\mathrm{data}^{-1}.
|
||||
\ee
|
||||
|
||||
Returning to the generic case where most energy-based models have \emph{intractable likelihoods}, we must estimate expectation values numerically. One way to do this is draw samples $\mathcal{S}_\mathrm{model}=\{ \mathbf{x}_i ' \}$ from the model $p_\theta(\mathbf{x})$ and evaluate expectation values using these samples:
|
||||
\be
|
||||
\langle h(\mathbf{x}) \rangle_\mathrm{model} = \int \mathrm{d}\mathbf{x} p_\theta(\mathbf{x}) h(\mathbf{x}) \approx \sum_{\mathbf{x}_i '\in \mathcal{S}_\mathrm{model}} g(\mathbf{x}_i').
|
||||
\ee
|
||||
The samples from the model $\mathbf{x}_i ' \in \mathcal{S}_\mathrm{model}$ are often referred to as \emph{fantasy particles} in the ML literature and can be generated using simple MCMC algorithms such as Metropolis-Hasting which are covered in most modern statistical physics classes. However, if the reader is unfamiliar with MCMC methods or wants a quick refresher, we recommend the concise and beautiful discussion of MCMC methods from both the physics and ML point-of-view in Chapters 29-32 of David MacKay's masterful book \cite{mackay2003information}.
|
||||
|
||||
|
||||
Finally, we note that once we have the fantasy particles (samples) from the model, we can also easily calculate the gradient of an arbitrary expectation value $\langle g(\mathbf{x}) \rangle_\mathrm{model}$ using what is commonly called the ``log-derivative trick'' in ML \cite{kleijnen1996optimization,fu2006gradient}:
|
||||
\bea
|
||||
{\partial \over \partial \theta_i}\langle g(\mathbf{x}) \rangle_\mathrm{model} &=& \int \mathrm{d}\mathbf{x} { \partial p_\theta(\mathbf{x}) \over \partial \theta_i} g(\mathbf{x}) \nonumber\\
|
||||
&=& \Big \langle {\partial \log{ p_\theta(\mathbf{x})} \over \partial \theta_i} g(\mathbf{x}) \Big \rangle_\mathrm{model} \nonumber \\
|
||||
&=& \langle O_i(\mathbf{x})g(\mathbf{x})\rangle_\mathrm{model}\nonumber \\
|
||||
&\approx& \sum_{\mathbf{x}_i' \in \mathcal{S}_\mathrm{model}} O_i(\mathbf{x})g(\mathbf{x}_i').
|
||||
\eea
|
||||
This expression allows us to take gradients of more complex cost functions beyond the MLE procedure discussed here.
|
||||
|
||||
|
||||
\subsection{Summary of the training procedure}
|
||||
|
||||
We now summarize the discussion above and present a general procedure for training an energy based model using SGD on the cost function (see Sec.~\ref{sec:gradient_descent}). Our goal is to fit the parameters of a model $p_{{\bf \lambda}}( \{\theta_i\} ) = Z^{-1} \mathrm e^{-E({\bf x}, \{\theta_i\})}$. Training the model involves the following steps:
|
||||
\begin{enumerate}
|
||||
\item Read a minibatch of data, $\{ {\bf x} \}$.
|
||||
\item Generate a random samplea (fantasy partscles) $\{ {\bf x}' \} \sim p_{{\bf \lambda}}$ using an MCMC algorithm (e.g., Metropolis-Hastings).
|
||||
\item Compute the gradient of log-likelihood using these samples and Eq.~\eqref{eq:Pos-neg-energy}, where the averages are taken over the minibatch of data and the fantasy particles/samples from the model, respectively.
|
||||
\item Use the gradient as input to one of the gradient based optimizers discussed in section Sec.~\ref{sec:gradient_descent}.
|
||||
\end{enumerate}
|
||||
In practice, it is helpful to supplement this basic procedure with some practical tricks that help training. As with discriminative neural networks, it is important to initialize the parameters properly and print summary statistics during the training procedure on the training and validation sets to prevent overfitting. These and many other little practical tricks have been nicely summarized in a short note from the Hinton group \cite{hinton2012practical}.
|
||||
|
||||
A major computational and practical limitation of these methods is that it is often hard to draw samples from generative models. MCMC methods often have long mixing-times (the time you have to run the Markov chain to get uncorrelated samples) and this can result in biased sampling. Luckily, we often do not need to know the gradients exactly for training ML models (recall that noisy gradient estimates often help the convergence of gradient descent algorithms) and we can significantly reduce the computational expense by running MCMC for a reasonable time window. We will exploit this observation extensively in the next section when we discuss how to train more complex energy-based models with hidden variables.
|
||||
|
||||
|
||||
\section{Deep Generative Models: Latent Variables and Restricted Boltzmann Machines (RBMs)}
|
||||
\label{sec:EnergyII}
|
||||
|
||||
The last section introduced many of the core ideas behind energy-based generative models. Here, we extend this discussion to energy-based models that include latent or hidden variables.
|
||||
|
||||
Including latent variables in generative models greatly enhances their expressive power -- allowing the model to represent sophisticated correlations between visible features without sacrificing trainability. By having multiple layers of latent variables, we can even construct powerful deep generative models that possess many of the same desirable properties as the deep, discriminative neural networks.
|
||||
|
||||
We begin with a discussion that tries to provide a simple intuition for why latent variables are such a powerful tool for generative models. Next, we introduce a powerful class of latent variable models called Restricted Boltzmann Machines (RBMs) and discuss techniques for training these models. After that, we discuss Deep Boltzmann Machines (DBMs), which have multiple layers of latent variables. We then introduce the new Paysage package for training energy-based models and demonstrate how to use it on the MNIST dataset and samples from the Ising model. We conclude by discussing recent physics literature related to energy-based generative models.
|
||||
|
||||
|
||||
\subsection{Why hidden (latent) variables?}
|
||||
|
||||
Latent or hidden variables are a powerful yet elegant way to encode sophisticated correlations between observable features. The underlying reason for this is that marginalizing over a subset of variables -- ``integrating out'' degrees of freedom in the language of physics -- induces complex interactions between the remaining variables. The idea that integrating out variables can lead to complex correlations is a familiar component of many physical theories. For example, when considering free electrons living on a lattice, integrating out phonons gives rise to higher-order electron-electron interactions (e.g.~superconducting or magnetic correlations). More generally, in the Wilsonian renormalization group paradigm, all effective field theories can be thought of as arising from integrating out high-energy degrees of freedom \cite{wilson1974renormalization}.
|
||||
|
||||
|
||||
Generative models with latent variables run this logic in reverse -- encode complex interactions between visible variables by introducing additional, hidden variables that interact with visible degrees of freedom in a simple manner, yet still reproduce the complex correlations between visible degrees in the data once marginalized over (integrated out). This allows us to encode complex higher-order interactions between the visible variables using simpler interactions at the cost of introducing new latent variables/degrees of freedom. This trick is also widely exploited in physics (e.g.~in the Hubbard-Stratonovich transformation \cite{hubbard1959calculation,stratonovich1957method} or the introduction of ghost fields in gauge theory \cite{faddeev1967feynman}).
|
||||
|
||||
To make these ideas more concrete, let us revisit the pairwise Ising model introduced in our discussion of MaxEnt models (see Eq.~\eqref{eq:pairwise-maxent}). The model is described
|
||||
by a Boltzmann distribution with energy
|
||||
\be
|
||||
E({\bf v}) = - \sum_i a_i v_i - \frac{1}{2} \sum_{ij} v_i J_{ij}v_j,
|
||||
\ee
|
||||
where $J_{ij}$ is a symmetric coupling matrix that encodes the pairwise constraints and $a_i$ enforce the single-variable constraint.
|
||||
|
||||
Our goal is to replace the complicated interactions between the visible variables $v_i$, encoded by $J_{ij}$, by interactions with a new set of latent variables $h_\mu$. In order to do this, it is helpful to rewrite the coupling matrix in a slightly different form. Using SVD, we can always express the coupling matrix in the form $J_{ij} = \sum_{\mu=1}^N W_{i \mu} W_{j \mu}$, where $\{W_{i\mu}\}_i$ are appropriately normalized singular vectors. In terms of $W_{i \mu}$, the energy takes the form
|
||||
\begin{equation}
|
||||
E_\mathrm{Hop}({\bf v}) = - \sum_i a_i v_i - \frac{1}{2} \sum_{ij\mu} v_i W_{i \mu} W_{j \mu} v_j.
|
||||
\label{eq:hopfield-energy}
|
||||
\end{equation}
|
||||
|
||||
We note that in the special case when both the $v_i \in \{-1, +1\}$ and $W_{i \mu} \in \{-1, +1 \}$ are binary variables, a model with this form of the energy function is known as the \emph{Hopfield model} \cite{hopfield1982neural, amit1985spin}. The Hopfield model has played an extremely important role in statistical physics, computational neuroscience, and machine learning, and a full discussion of its properties is well beyond the scope of this review (see \cite{amit1992modeling} for a beautiful discussion that combines all these perspectives). Therefore, here we refer to all energy functions of the form Eq.~\eqref{eq:hopfield-energy} as (generalized) Hopfield models, even for the case when the $W_{i \mu}$ are continuous variables.
|
||||
|
||||
We now ``decouple'' the $v_i$ by introducing a set of normally, distributed continuous latent variables $h_{\mu}$ (in condensed matter physics this is called a Hubbard-Stratonovich transformation). Using the usual identity for Gaussian integrals, we can rewrite the Boltzmann distribution for the generalized Hopfield model as
|
||||
\bea
|
||||
p(\mathbf{v}) &=& {\mathrm e^{ \sum_i a_i v_i + \frac{1}{2} \sum_{ij\mu} v_i W_{i \mu} W_{j \mu} v_j} \over Z} \nonumber \\
|
||||
&=& { \mathrm e^{\sum_{i} a_i v_i }\prod_\mu \int \mathrm{d}h_\mu \mathrm e^{- \frac{1}{2} \sum_{\mu} h_{\mu}^2 - \sum_{i} v_i W_{i \mu} h_{\mu}} \over Z} \nonumber \\
|
||||
&= &\int \mathrm{d}\mathbf{h} \mathrm e^{-E({\mathbf v}, {\mathbf h})} \over Z
|
||||
\label{eq:HS-hopfield}
|
||||
\eea
|
||||
where $E( {\bf v}, {\bf h} )$ is a joint energy functional of both the latent and visible variables of the form
|
||||
\be
|
||||
E( {\bf v}, {\bf h} ) = - \sum_{i} a_i v_i + \frac{1}{2} \sum_{\mu} h_{\mu}^2 - \sum_{i \mu} v_i W_{i \mu} h_{\mu}.
|
||||
\label{eq:hopfield-rbm}
|
||||
\ee
|
||||
We can also use the energy function $E({\mathbf v}, {\mathbf h})$ to define a new energy-based model $p({\bf v}, {\bf h})$ on both the latent and visible variables
|
||||
\be
|
||||
p({\bf v}, {\bf h})={\mathrm e^{-E({\mathbf v}, {\mathbf h})} \over Z'}.
|
||||
\ee
|
||||
Marginalizing over latent variables of course gives us back the generalized Hopfield model \cite{barra2012equivalence}
|
||||
\be
|
||||
p(\mathbf{v})= \int d{\bf h} p({\bf v}, {\bf h})={\mathrm e^{-E_\mathrm{Hop}({\bf v})} \over Z}.
|
||||
\ee
|
||||
|
||||
|
||||
Notice that $E( {\bf v}, {\bf h} )$ contains no direct interactions between visible degrees of freedom (or between hidden degree of freedom). Instead, the complex correlations between the $v_i$ are encoded in the interaction between the visible $v_i$ and latent variables $h_\mu$. It turns out that the model presented here is a special case of a more general class of powerful energy-based models called Restricted Boltzmann Machines (RBMs).
|
||||
|
||||
|
||||
\subsection{Restricted Boltzmann Machines (RBMs)}
|
||||
\label{sec:RBMs}
|
||||
|
||||
\begin{figure}[t]
|
||||
\includegraphics[width=1.0\columnwidth]{Energy/RBM.pdf}
|
||||
\caption{A Restricted Boltzmann Machine (RBM) consists of visible units $v_i$ and hidden units $h_\mu$ that interact
|
||||
with each other through interactions of the form $W_{i\mu}v_i h_\mu$. Importantly, there are no interactions between visible units themselves or hidden units themselves. }
|
||||
\label{fig:RBM-architecture}
|
||||
\end{figure}
|
||||
|
||||
A Restricted Boltzmann Machine (RBM) is an energy-based model with both visible and hidden units where the visible and hidden units interact with each other but do not interact among themselves. The energy function of an RBM takes the general functional form
|
||||
\begin{equation}
|
||||
E({\bf v}, {\bf h}) = -\sum_i a_i( v_i) - \sum_{\mu} b_{\mu}(h_{\mu}) - \sum_{i \mu} W_{i \mu} v_i h_{\mu}
|
||||
\label{eq:rbm}
|
||||
\end{equation}
|
||||
where $a_i(\cdot)$ and $b_{\mu}(\cdot)$ are functions that we are free to choose. The most common choice is:
|
||||
\[
|
||||
a_i( v_i):=
|
||||
\begin{cases}
|
||||
a_i v_i,& \text{if } v_i \in \{ 0, 1\} \text{ is binary} \\
|
||||
{v_i^2 \over 2 \sigma_i^2}, & \text{if } v_i \in \mathbb{R} \text{ is continuous},
|
||||
\end{cases}
|
||||
\]
|
||||
and
|
||||
\[
|
||||
b_\mu( h_\mu):=
|
||||
\begin{cases}
|
||||
b_\mu h_\mu ,& \text{if } h_\mu \in \{ 0, 1\} \text{ is binary} \\
|
||||
{h_\mu^2 \over 2 \sigma_\mu^2}, & \text{if } h_\mu \in \mathbb{R} \text{ is continuous}.
|
||||
\end{cases}
|
||||
\]
|
||||
For this choice of $a_i(\cdot)$ and $b_{\mu}(\cdot)$, layers consisting of discrete binary units are often called Bernoulli layers and layers consisting of continuous variables
|
||||
are often called Gaussian layers. The basic bipartite structure of an RBM -- i.e., a visible and hidden layer that interact with each other but not among themselves -- is often depicted using a graph of the form shown in Fig. \ref{fig:RBM-architecture}.
|
||||
|
||||
|
||||
An RBM can have different properties depending on whether the hidden and visible layers are taken to be Bernoulli or Gaussian. The most common choice is to have both the visible and hidden units be Bernoulli. This is what is typically meant by an RBM. However, other combinations are also possible and used in the ML literature. When all the units are continuous, the RBM reduces to a multi-dimensional Gaussian with a very particular correlation structure. When the hidden units are continuous and the visible units are discrete, the RBM is equivalent to a generalized Hopfield model (see discussion above). When the the visible units are continuous and the hidden units are discrete, the RBM is often called a Gaussian Bernoulli Restricted Boltzmann Machine \cite{dahl2010phone, hinton2006reducing}. It is even possible to perform multi-modal learning with a mixture of continuous and discrete variables. For all these architectures, the important point is that all interactions occur only between the visible and hidden units and there are no interactions between units within the hidden or visible layers (see Fig. ~\ref{fig:RBM-architecture}). This is analogous to Quantum Electrodynamics, where a free fermion and a free photon interact with one another but not among themselves.
|
||||
|
||||
Specifying a generative model with this bipartite interaction structure has two major advantages: (i) it enables capturing both pairwise {\it and higher-order} correlations between the visible units and (ii) it makes it easier to sample from the model using an MCMC method known as block Gibbs sampling, which in turn makes the model easier to train.
|
||||
|
||||
Before discussing training, it is worth better understanding the kind of correlations that can be captured using an RBM. To do so, we can marginalize over the hidden units and ask about the resulting distribution over just the visible units
|
||||
\be
|
||||
p({\bf v})=\int \mathrm{d}{\bf h} p(\mathbf{v}, \mathbf{h}) =\int \mathrm{d}{\bf h} {e^{-E({\bf v}, {\bf h})} \over Z}
|
||||
\ee
|
||||
(where the integral should be replaced by a trace in all expressions for discrete units).
|
||||
|
||||
We can also define a marginal energy using the expression
|
||||
\be
|
||||
p({\bf v}):={\mathrm e^{-E({\bf v})} \over Z}.
|
||||
\ee
|
||||
Combining these equations,
|
||||
\begin{align}
|
||||
&E({\bf v}) = - \log \int \mathrm d {\bf h} \mathrm e^{-E({\bf v}, {\bf h})} \nonumber \\
|
||||
&= - \sum_i a_i (v_i) - \sum_{\mu} \log \int \mathrm d h_{\mu} \mathrm e^{ b_{\mu} (h_{\mu}) + \sum_i v_i W_{i \mu} h_{\mu}} \nonumber
|
||||
\end{align}
|
||||
To understand what correlations are captured by $p({\bf v})$ it is helpful to introduce the distribution
|
||||
\be
|
||||
q_{\mu}(h_{\mu}) = { \mathrm e^{ b_{\mu}(h_{\mu})} \over Z}
|
||||
\ee
|
||||
of hidden units $h_{\mu}$, ignoring the interactions between ${\bf v}$ and ${\bf h}$, and the cumulant generating function
|
||||
\begin{equation}
|
||||
K_{\mu}(t) := \log \int \mathrm d h_{\mu} q_{\mu}(h_{\mu}) \mathrm e^{t h_{\mu}} = \sum_{n} \kappa_{\mu}^{(n)} \frac{t^n}{n!}.
|
||||
\end{equation}
|
||||
$K_{\mu}(t)$ is defined such that the $n^{th}$ cumulant is $\kappa_{\mu}^{(n)} = \partial_t^{n} K_{\mu}\vert_{t=0}$.
|
||||
|
||||
The cumulant generating function appears in the marginal free-energy of the visible units, which can be rewritten (up to a constant term) as:
|
||||
\begin{align}
|
||||
&E({\bf v}) = - \sum_i a_i (v_i) - \sum_{\mu} K_{\mu}\left(\sum_i W_{i \mu} v_i\right) \nonumber \\
|
||||
&= - \sum_i a_i (v_i) - \sum_{\mu} \sum_{n} \kappa_{\mu}^{(n)} \frac{(\sum_i W_{i \mu} v_i)^n}{n!} \nonumber \\
|
||||
&= - \sum_i a_i (v_i) - \sum_i \left( \sum_{\mu} \kappa_{\mu}^{(1)} W_{i \mu} \right) v_i \nonumber \\
|
||||
&- \frac{1}{2} \sum_{ij} \left(\sum_{\mu} \kappa_{\mu}^{(2)}W_{i \mu} W_{j \mu} \right) v_i v_j + \ldots
|
||||
\label{eq:cumulant-rbm}
|
||||
\end{align}
|
||||
We see that the marginal energy includes all orders of interactions between the visible units, with the $n$-th order cumulants of $q_{\mu}(h_{\mu})$ weighting the $n$-th order interactions between the visible units. In the case of the Hopfield model we discussed previously, $q_{\mu}(h_{\mu})$ is a standard Gaussian distribution where the mean is $\kappa_{\mu}^{(1)} = 0$, the variance is $\kappa_{\mu}^{(2)} = 1$, and all higher order cumulants are zero. Plugging these cumulants into Eq.~\eqref{eq:cumulant-rbm} recovers Eq.~\eqref{eq:hopfield-rbm}.
|
||||
|
||||
These calculations make clear the underlying reason for the incredible representational power of RBMs with a Bernoulli hidden layer. Each hidden unit can encode interactions of arbitrarily high order. By combining many different hidden units, we can encode very complex interactions at all orders. Moreover, we can learn which order of correlations/interactions are important directly from the data instead of having to specify them ahead of time as we did in the MaxEnt models. This highlights the power of generative models with even the simplest interactions between visible and latent variables to encode, learn, and represent complex correlations present in the data.
|
||||
|
||||
\subsection{Training RBMs}
|
||||
|
||||
|
||||
RBMs are a special class of energy-based generative models, which can be trained using the Maximum Likelihood Estimation (MLE) procedure described in detail in Sec.~\ref{sec:EnergyI}. To briefly recap, first, we must choose a cost function -- for MLE this is just the negative log-likelihood with or without an additional regularization term to prevent overfitting. We then minimize this cost function using one of the Stochastic Gradient Descent (SGD) methods described in Sec.~\ref{sec:gradient_descent}.
|
||||
|
||||
|
||||
The gradient itself can be calculated using Eq.~\eqref{eq:Pos-neg-energy}. For example, for the Bernoulli-Bernoulli RBM in Eq.~\eqref{eq:rbm} we have
|
||||
\bea
|
||||
{\partial \mathcal{L}(\{W_{i \mu}, a_i, b_\mu \}) \over \partial W_{i \mu}} &=& \langle v_i h_{\mu} \rangle_\mathrm{data} - \langle v_i h_{\mu} \rangle_\mathrm{model} \\
|
||||
{\partial \mathcal{L}(\{W_{i \mu}, a_i, b_\mu \}) \over \partial a_i} &=& \langle v_i \rangle_\mathrm{data} - \langle v_i \rangle_\mathrm{model} \nonumber \\
|
||||
{\partial \mathcal{L}(\{W_{i \mu}, a_i, b_\mu \}) \over \partial b_\mu} &=& \langle h_\mu \rangle_\mathrm{data} - \langle h_\mu \rangle_\mathrm{model},\nonumber
|
||||
\eea
|
||||
where the positive expectation with respect to data is understood to mean sampling from the model while clamping the visible units to their observed values in the data. As before, calculating the negative phase of the gradient (the expectation value with respect to the model) requires that we draw samples from the model. Luckily, the bipartite form of the interactions in RBMs were specifically chosen with this in mind.
|
||||
|
||||
\subsubsection{Gibbs sampling and contrastive divergence (CD)}
|
||||
|
||||
\begin{figure}[t]
|
||||
\includegraphics[width=1.0\columnwidth]{Energy/ContrastiveDivergence.pdf}
|
||||
\caption{ (Top) To draw fantasy particles (samples from model) we can perform alternating (block) Gibbs sampling between the visible and hidden layers starting
|
||||
with a sample from the data using the marginal distributions $p(\mathbf{h}|\mathbf{v})$ and $p(\mathbf{v}| \mathbf{h})$. The ``time'' $t$ corresponds to the time in the
|
||||
Markov chain for the Monte Carlo and measures the number of passes between the visible and hidden states. (Middle) In Contrastive Divergence (CD), we approximately
|
||||
sample the model by terminating the Gibbs sampling after $n$ steps (CD-n) starting from the data. (C) In Persistent Contrastive Divergence (PCD), instead of restarting
|
||||
our sample from the data, we initialize the sampler with the fantasy particles calculated from the model at the last SGD step.}
|
||||
\label{fig:CD-PCD}
|
||||
\end{figure}
|
||||
|
||||
|
||||
|
||||
|
||||
The bipartite interaction structure of an RBM makes it possible to calculate expectation values using a Markov Chain Monte Carlo (MCMC) method known as Gibbs sampling. The key reason for this is that since there are no interactions of visible units with themselves or hidden units with themselves; the visible and hidden units of an RBM are conditionally independent:
|
||||
\begin{align}
|
||||
p({\bf v} | {\bf h}) &= \prod_i p(v_i |{\bf h}) \nonumber \\
|
||||
p({\bf h} | {\bf v}) &= \prod_\mu p(h_\mu|{\bf v}),
|
||||
\end{align}
|
||||
with
|
||||
\begin{eqnarray}
|
||||
\label{eq:RBM-conditional}
|
||||
p(v_i = 1 | {\bf h}) &=& \sigma(a_i + \sum_{\mu} W_{i \mu} h_{\mu}) \\
|
||||
p(h_{\mu} = 1 | {\bf v}) &=& \sigma(b_{\mu}+ \sum_{i} W_{i \mu} v_{i}) \nonumber
|
||||
\end{eqnarray}
|
||||
and where $\sigma(x) = 1 / (1 + \mathrm e^{-x})$ is the sigmoid function.
|
||||
|
||||
Using these expressions it is easy to compute expectation values with respect to the data. The input to gradient descent is a minibatch of observed data. For each sample in the minibatch, we simply clamp the visible units to the observed values and apply Eq.~\eqref{eq:RBM-conditional} using the probability for the hidden variables. We then average over all samples in the minibatch to calculate expectation values with respect to the data. To calculate expectation values with respect to the model, we use (block) Gibbs sampling. The idea behind (block) Gibbs sampling is to iteratively sample from the conditional distributions ${\bf h}_{t+1} \sim p({\bf h} | {\bf v}_{t})$ and ${\bf v}_{t+1} \sim p({\bf v} | {\bf h}_{t+1})$ (see Figure \ref{fig:CD-PCD}, top). Since the units are conditionally independent, each step of this iteration can be performed by simply drawing random numbers. The samples are guaranteed to converge to the equilibrium distribution of the model in the limit that $t \rightarrow \infty$. At the end of the Gibbs sampling procedure, one ends up with a minibatch of samples (fantasy particles).
|
||||
|
||||
One drawback of Gibbs sampling is that it may take many back and forth iterations to draw an independent sample. For this reason, the Hinton group introduced an approximate Gibbs sampling technique called Contrastive Divergence (CD) \cite{hinton2002training, hinton2006fast}. In CD-$n$, we just perform $n$ iterations of (block) Gibbs sampling, with $n$ often taken to be as small as 1 (see Figure \ref{fig:CD-PCD})! The price for this truncation is, of course, that we are not drawing samples from the true model distribution. But for our purpose -- using the expectations to estimate the gradient for SGD -- CD-$n$ has been proven to work reasonably well. As long as the approximate gradients are reasonably correlated with the true gradient, SGD will move in a reasonable direction. CD-$n$ of course does come at a price. Truncating the Gibbs sampler prevents sampling far away from the starting point, which for CD-$n$ are the data points in the minibatch. Therefore, our generative model will be much more accurate around regions of feature space close to our training data. Thus, as is often the case in ML, CD-$n$ sacrifices the ability to generalize to some extent in order to make the model easier to train.
|
||||
|
||||
Some of these undesirable features can be tempered by using a slightly different variant of CD called Persistent Contrastive Divergence (PCD) \cite{tieleman2009using}. In PCD, rather than restarting the Gibbs sampler from the data at each gradient descent step, we start the Gibbs sampling at the fantasy particles (samples from the model) in the last gradient descent step (see Figure \ref{fig:CD-PCD}). Since parameters change slowly compared to the Gibbs sampling, samples that are high probability at one step of the SGD are also likely to be high probability at the next step. This ensures that PCD does not introduce large errors in the estimation of the gradients. The advantage of using fantasy particles to initialize the Gibbs sampler is to allow PCD to explore parts of the feature space that are much further from the training dataset than once could reach with ordinary CD.
|
||||
|
||||
|
||||
|
||||
\subsection{Deep Boltzmann Machine}
|
||||
|
||||
In this section, we introduce Deep Boltzmann Machines (DBMs). Unlike RBMs, DBMs possess multiple hidden layers and were the first models rebranded as ``deep learning'' \cite{hinton2006fast,hinton2006reducing} \footnote{Technically, these were Deep Belief Networks (DBNs) where only the top layer was undirected}. Many of the advantages that are thought to stem from having deep layers were already discussed in Sec.~\ref{sec:DNN-III} in the context of discriminative DNNs. Here, we revisit many of the same themes with emphasis on energy-based models.
|
||||
|
||||
An RBM is composed of two layers of neurons that are connected via an undirected graph, see Fig.~\ref{fig:RBM-architecture}. As a result, it is possible to perform sampling ${\bf v} \sim p({\bf v} | {\bf h})$ and inference ${\bf h} \sim p({\bf h} |{\bf v})$ with the same model. As with the Hopfield model, we can view each of the hidden units as representative of a pattern, or feature, that could be present in the data. (In general, one should think of activity patterns of hidden units representing features in the data.) The inference step involves assigning a probability to each of these features that expresses the degree to which each feature is present in a given data sample. In an RBM, hidden units do not influence each other during the inference step, i.e.~hidden units are conditionally independent given the visible units. There are a number of reasons why this is unsatisfactory. One reason is the desire for sparse, distributed representations, where each observed visible vector will strongly activate a few (i.e.~more than one but only a very small fraction) of the hidden units. In the brain, this is thought to be achieved by inhibitory lateral connections between neurons. However, adding lateral intra-layer connections between the hidden units makes the distribution difficult to sample from, so we need to come up with another way of creating connections between the hidden units.
|
||||
|
||||
With the Hopfield model, we saw that pairwise linear connections between neurons can be mediated through another layer. Therefore, a simple way to allow for effective connections between the hidden units is to add another layer of hidden units. Rather than just having two layers, one visible and one hidden, we can add additional layers of latent variables to account for the correlations between hidden units. Ideally, as one adds more and more layers, one might hope that the correlations between hidden variables become smaller and smaller deeper into the network. This basic logic is reminiscent of renormalization procedures that seek to decorrelate layers at each step \cite{mehta2014exact, vidal2007entanglement, li2018neural}. The price of adding additional layers is that the models become harder to train.
|
||||
|
||||
|
||||
Training DBMs is more subtle than RBMs due to the difficulty of propagating information from visible to hidden units. However, Hinton and collaborators realized that some of these problems could be alleviated via a layerwise procedure. Rather than attempting to the train the whole DBM at once, we can think of the DBM as a stack of RBMs (see Fig.~\ref{fig:DBM}). One first trains the bottom two layers of the DBM -- treating it as if it is a stand-alone RBM. Once this bottom RBM is trained, we can generate ``samples'' from the hidden layer and use these samples as an input to the next RBM (consisting of the first and second hidden layer -- purple hexagons and green squares in Fig. \ref{fig:DBM}). This procedure can then be repeated to pretrain all layers of the DBM.
|
||||
|
||||
|
||||
This pretraining initializes the weights so that SGD can be used effectively when the network is trained in a supervised fashion. In particular, the pretraining helps the gradients to stay well behaved rather than vanish or blow up -- a problem that we discussed extensively in the earlier sections of DNNs. It is worth noting that once pretrained, we can use the usual Boltzmann learning rules (Eq.~\eqref{eq:Pos-neg-energy}) to fine-tune the weights and improve the performance of the DBM \cite{ hinton2006fast,hinton2006reducing}. As we demonstrate in the next section, the Paysage package presented here can be used to both construct and train DBNs using such a pretraining procedure.
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -4,23 +4,114 @@ DATE: today
|
||||
|
||||
|
||||
!split
|
||||
===== What is Machine Learning? =====
|
||||
===== Introduction =====
|
||||
|
||||
Add definition of definition of machine learning
|
||||
|
||||
Machine learning is the science of giving computers the ability to
|
||||
learn without being explicitly programmed. The idea is that there
|
||||
exist generic algorithms which can be used to find patterns in a broad
|
||||
class of data sets without having to write code specifically for each
|
||||
problem. The algorithm will build its own logic based on the data.
|
||||
Statistics, data science and machine learning form important fields of
|
||||
research in modern science. They describe how to learn and make
|
||||
predictions from data, as well allowing us to extract important
|
||||
correlations about physical process and the underlying laws of motion
|
||||
in large data sets. The latter, big data sets, appear
|
||||
frequently in essentially all disciplines, from the traditional Science,
|
||||
Technology, Mathematics and Engineering fields to Life Science, Law, education research,
|
||||
the Humanities and
|
||||
the Social Sciences. It has become more and more common to see
|
||||
research projects on big data in for example the Social
|
||||
Sciences where extracting patterns from complicated survey data is one of many research directions.
|
||||
Having a solid grasp of data analysis and machine learning
|
||||
is thus becoming central to scientific computing in many
|
||||
fields, and competences and skills within the fields of machine learning
|
||||
and scientific computing are nowadays strongly requested by many
|
||||
potential employers. The latter cannot be overstated, familiarity with
|
||||
machine learning has almost become a prerequisite for many of the most
|
||||
exciting employment opportunities, whether they are in bioinformatics,
|
||||
life science, physics or finance, in the private or the public
|
||||
sector. This author has had several students or met students who have
|
||||
been hired recently based on their skills and competences in
|
||||
scientific computing and data science, often with marginal knowledge
|
||||
of machine learning.
|
||||
|
||||
Machine learning is a subfield of computer science, and is closely
|
||||
related to computational statistics. It evolved from the study of
|
||||
pattern recognition in artificial intelligence (AI) research, and has
|
||||
made contributions to AI tasks like computer vision, natural language
|
||||
processing and speech recognition. It has also, especially in later
|
||||
years, found applications in a wide variety of other areas, including
|
||||
bioinformatics, economy, physics, finance and marketing.
|
||||
processing and speech recognition.
|
||||
Machine learning represents the
|
||||
science of giving computers the ability to learn without being
|
||||
explicitly programmed. The idea is that there exist generic
|
||||
algorithms which can be used to find patterns in a broad class of data
|
||||
sets without having to write code specifically for each problem. The
|
||||
algorithm will build its own logic based on the data.
|
||||
|
||||
Machine learning is an extremely rich field, in spite of its young age. The
|
||||
increases we have seen during the last three decades in computational
|
||||
capabilities have been followed by developments of methods and
|
||||
techniques for analyzing and handling large date sets, relying heavily
|
||||
on statistics, computer science and mathematics. The field is rather
|
||||
new and developing rapidly. Popular software packages written in
|
||||
Python for machine learning like "Scikit-learn":"http://scikit-learn.org/stable/", "Tensorflow":"https://www.tensorflow.org/",
|
||||
"PyTorch":"http://pytorch.org/" and "Keras":"https://keras.io/", all freely available at their respective GitHub sites,
|
||||
encompass communities of developers in the thousands or more. And the number
|
||||
of code developers and contributors keeps increasing. Not all the
|
||||
algorithms and methods can be given a rigorous mathematical
|
||||
justification, opening up thereby large rooms for experimenting
|
||||
and trial and error and thereby exciting new developments.
|
||||
However, a solid command of linear algebra, multivariate theory,
|
||||
probability theory, statistical data analysis,
|
||||
understanding errors and Monte Carlo methods are central elements in a proper understanding of many of
|
||||
algorithms and methods we will discuss.
|
||||
|
||||
!split
|
||||
===== Learning outcomes =====
|
||||
|
||||
These lectures aim at giving you an overview of central aspects of
|
||||
statistical data analysis as well as some of the central algorithms
|
||||
used in machine learning. We will introduce a variety of central
|
||||
algorithms and methods essential for studies of data analysis and
|
||||
machine learning.
|
||||
|
||||
Hands-on projects and experimenting with data and algorithms plays a central role in
|
||||
these lectures, and our hope is, through the various
|
||||
projects and exercies, to expose you to fundamental
|
||||
research problems in these fields, with the aim to reproduce state of
|
||||
the art scientific results. You will learn to develop and
|
||||
structure large codes for studying these systems, get acquainted with
|
||||
computing facilities and learn to handle large scientific projects. A
|
||||
good scientific and ethical conduct is emphasized throughout the
|
||||
course. More specifically, you will
|
||||
|
||||
o learn about basic data analysis, Bayesian statistics, Monte Carlo methods, data optimization and machine learning;
|
||||
o be capable of extending the acquired knowledge to other systems and cases;
|
||||
o Have an understanding of central algorithms used in data analysis and machine learning;
|
||||
o Gain knowledge of central aspects of Monte Carlo methods, Markov chains, Gibbs samplers and their possible applications, from numerical integration to simulation of stock markets;
|
||||
o Understand methods for regression and classification;
|
||||
o Learn about neural network, genetic algorithms and Boltzmann machines;
|
||||
o Work on numerical projects to illustrate the theory. The projects play a central role and you are expected to know modern programming languages like Python or C++, in addition to a basic knowledge of linear algebra (typically taught during the first one or two years of undergraduate studies).
|
||||
|
||||
There are several topics we will cover here, spanning from a
|
||||
statistical data analysis and its basic concepts such expectation
|
||||
values, variance, covariance, correlation functions and errors, via
|
||||
well-known probability distribution functions like uniform
|
||||
distribution, the binomial distribution, the Poisson distribution and
|
||||
simple and multivariate normal distributions to central elements of
|
||||
Bayesian statistics and modeling. We will also remind the reader about
|
||||
central elements from linear algebra and standard methods based on
|
||||
linear algebra used to fit functions such Cubic splines and gradient
|
||||
methods for data optimization and the Singular-value decomposition and
|
||||
least square methods for parameterizing data.
|
||||
|
||||
We will also cover Monte Carlo methods, Markov chains, well-known
|
||||
algorithms for sampling stochastic events like the Metropolis-Hastings
|
||||
and Gibbs sampling methods. An important aspect of all our
|
||||
calculations is a proper estimation of errors. Here we will also
|
||||
discuss famous resampling techniques like the blocking, bootstrapping
|
||||
and jackknife methods.
|
||||
|
||||
The second part of the material covers several algorithms used in
|
||||
machine learning.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Types of Machine Learning =====
|
||||
@@ -31,7 +122,7 @@ In *supervised learning* we know the answer to a problem,
|
||||
and let the computer deduce the logic behind it. On the other hand, *unsupervised learning*
|
||||
is a method for finding patterns and relationship in data sets without any prior knowledge of the system.
|
||||
Some authours also operate with a third category, namely *reinforcement learning*. This is a paradigm
|
||||
of learning inspired by behavioural psychology, where learning is achieved by trial-and-error,
|
||||
of learning inspired by behavioral psychology, where learning is achieved by trial-and-error,
|
||||
solely from rewards and punishment.
|
||||
|
||||
Another way to categorize machine learning tasks is to consider the desired output of a system.
|
||||
@@ -44,27 +135,52 @@ Some of the most common tasks are:
|
||||
* Clustering: Data are divided into groups with certain common traits, without knowing the different groups beforehand. It is thus a form of unsupervised learning.
|
||||
|
||||
|
||||
!split
|
||||
===== Different algorithms =====
|
||||
In this course we will build our machine learning approach on a statistical foundation, with elements
|
||||
from data analysis, stochastic processes etc before we proceed with the following machine learning algorithms
|
||||
The methods we cover have three main topics in common, irrespective of
|
||||
whether we deal with supervised or unsupervised learning. The first
|
||||
ingredient is normally our data set, the second is a model which is
|
||||
normally a function of some parameters. The last ingredient is a
|
||||
so-called _cost_ function which allows us to present an estimate on
|
||||
how good our model is in reproducing the data it is supposed to train.
|
||||
|
||||
o Linear regression and its variants
|
||||
Here we will build our machine learning approach on elements of the
|
||||
statistical foundation discussed above, with elements from data
|
||||
analysis, stochastic processes etc. We will discuss the following
|
||||
machine learning algorithms
|
||||
|
||||
o Linear regression and its variants, in essence polynomial regression
|
||||
o Decision tree algorithms, from simpler to more complex ones
|
||||
o Nearest neighbors models
|
||||
o Bayesian statistics
|
||||
o Bayesian statistics and regression
|
||||
o Support vector machines and finally various variants of
|
||||
o Artifical neural networks
|
||||
o Artifical neural networks and deep learning
|
||||
|
||||
Before we proceed however, there are several practicalities with data analysis and software tools we would
|
||||
like to present. These tools will help us in our understanding of various machine learning algorithms.
|
||||
|
||||
Our emphasis here is on understanding the mathematical aspects of different algorithms, however, where possible
|
||||
we will emphasize the importance of using available software.
|
||||
Before we proceed however, there are several practicalities with data
|
||||
analysis and software tools we would like to present. These tools will
|
||||
help us in our understanding of various machine learning algorithms.
|
||||
|
||||
Our emphasis here is on understanding the mathematical aspects of
|
||||
different algorithms, however, where possible we will emphasize the
|
||||
importance of using available software. We start thus with a hands-on
|
||||
and top-down approach machine learning. The aim is thus to start with
|
||||
relevant data and use these to introduce statistical data analysis
|
||||
concepts and machine learning algorithms before we delve into the
|
||||
algorithms themselves. The examples we will use start with a simple
|
||||
third-order polynomial with random noise added, and using the Python
|
||||
software package "Scikit-learn":"http://scikit-learn.org/stable/" we
|
||||
will introduce various machine learning algorithm s to make fits of
|
||||
the data data and predictions. We move thereafter to more interesting
|
||||
cases such as the simulation of financial transactions or disease
|
||||
models. These are examples where we can easily set up the data and
|
||||
then use machine learning algorithms using included in for example _scikit-learn_. Another model we
|
||||
will consider is the so-called Ising model. Here we will use this
|
||||
model to produce data for selected spin configurations and attempt to classify the data.
|
||||
Finally, our last example consists of economic data from the OECD.
|
||||
|
||||
|
||||
!split
|
||||
===== Software and needed installations =====
|
||||
|
||||
We will make intensive use of python as programming language and the myriad of available libraries.
|
||||
Furthermore, you will find IPython/Jupyter notebooks invaluable in your work.
|
||||
You can run _R_ codes in the Jupyter/IPython notebooks, with the immediate benefit of visualizing your data.
|
||||
@@ -107,7 +223,7 @@ To install _R_ with Jupyter notebook "following the link here":"https://mpacer.o
|
||||
===== Installing R, C++, cython or Julia =====
|
||||
|
||||
|
||||
For the C++ affecianodas, Jupyter/IPython notebook allows you also to install C++ and run codes written in this language
|
||||
For the C++ aficionados, Jupyter/IPython notebook allows you also to install C++ and run codes written in this language
|
||||
interactively in the browser. Since we will emphasize writing many of the algorithms yourself, you can thus opt for
|
||||
either Python or C++ as programming languages.
|
||||
|
||||
@@ -166,11 +282,145 @@ line = np.linspace(-3,3,1000,endpoint=False).reshape(-1,1)
|
||||
reg = DecisionTreeRegressor(min_samples_split=3).fit(x,y)
|
||||
plt.plot(line, reg.predict(line), label="decision tree")
|
||||
regline = LinearRegression().fit(x,y)
|
||||
plt.plot(line, regline.predict(line), label= "Linear Rgression")
|
||||
plt.plot(line, regline.predict(line), label= "Linear Regression")
|
||||
plt.show()
|
||||
!ec
|
||||
!eblock
|
||||
|
||||
|
||||
!split
|
||||
===== Simple regression model =====
|
||||
Add info about the equations
|
||||
!bc pycod
|
||||
# Importing various packages
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
x = 2*np.random.rand(100,1)
|
||||
y = 4+3*x+np.random.randn(100,1)
|
||||
|
||||
xb = np.c_[np.ones((100,1)), x]
|
||||
theta = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
|
||||
xnew = np.array([[0],[2]])
|
||||
xbnew = np.c_[np.ones((2,1)), xnew]
|
||||
ypredict = xbnew.dot(theta)
|
||||
|
||||
plt.plot(xnew, ypredict, "r-")
|
||||
plt.plot(x, y ,'ro')
|
||||
plt.axis([0,2.0,0, 15.0])
|
||||
plt.xlabel(r'$x$')
|
||||
plt.ylabel(r'$y$')
|
||||
plt.title(r'Linear Regression')
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Simple regression model, now using _scikit-learn_ =====
|
||||
Add info about the equations
|
||||
!bc pycod
|
||||
# Importing various packages
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.linear_model import LinearRegression
|
||||
|
||||
x = 2*np.random.rand(100,1)
|
||||
y = 4+3*x+np.random.randn(100,1)
|
||||
linreg = LinearRegression()
|
||||
linreg.fit(x,y)
|
||||
xnew = np.array([[0],[2]])
|
||||
ypredict = linreg.predict(xnew)
|
||||
|
||||
plt.plot(xnew, ypredict, "r-")
|
||||
plt.plot(x, y ,'ro')
|
||||
plt.axis([0,2.0,0, 15.0])
|
||||
plt.xlabel(r'$x$')
|
||||
plt.ylabel(r'$y$')
|
||||
plt.title(r'Random numbers ')
|
||||
plt.show()
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Simple regression model with gradient descent =====
|
||||
Add info about the equations, play around with different learning rates
|
||||
!bc pycod
|
||||
# Importing various packages
|
||||
from math import exp, sqrt
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
x = 2*np.random.rand(100,1)
|
||||
y = 4+3*x+np.random.randn(100,1)
|
||||
|
||||
xb = np.c_[np.ones((100,1)), x]
|
||||
theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
|
||||
print(theta_linreg)
|
||||
theta = np.random.randn(2,1)
|
||||
|
||||
eta = 0.1
|
||||
Niterations = 1000
|
||||
m = 100
|
||||
|
||||
for iter in range(Niterations):
|
||||
gradients = 2.0/m*xb.T.dot(xb.dot(theta)-y)
|
||||
theta -= eta*gradients
|
||||
|
||||
print(theta)
|
||||
xnew = np.array([[0],[2]])
|
||||
xbnew = np.c_[np.ones((2,1)), xnew]
|
||||
ypredict = xbnew.dot(theta)
|
||||
ypredict2 = xbnew.dot(theta_linreg)
|
||||
plt.plot(xnew, ypredict, "r-")
|
||||
plt.plot(xnew, ypredict2, "b-")
|
||||
plt.plot(x, y ,'ro')
|
||||
plt.axis([0,2.0,0, 15.0])
|
||||
plt.xlabel(r'$x$')
|
||||
plt.ylabel(r'$y$')
|
||||
plt.title(r'Random numbers ')
|
||||
plt.show()
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Simple regression model with stochastic gradient descent =====
|
||||
Add info about the equations, play around with different learning rates
|
||||
!bc pycod
|
||||
# Importing various packages
|
||||
from math import exp, sqrt
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.linear_model import SGDRegressor
|
||||
|
||||
x = 2*np.random.rand(100,1)
|
||||
y = 4+3*x+np.random.randn(100,1)
|
||||
|
||||
xb = np.c_[np.ones((100,1)), x]
|
||||
theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
|
||||
print(theta_linreg)
|
||||
sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
|
||||
sgdreg.fit(x,y.ravel())
|
||||
print(sgdreg.intercept_, sgdreg.coef_)
|
||||
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== Polynomial regression =====
|
||||
!bc pycod
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Predator-Prey model from ecology =====
|
||||
|
||||
@@ -500,4 +750,294 @@ text(0, 0.5, paste("y =x^ (", power, " +/- ", power.se, ")", sep = ""), pos = 4)
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Example: ecoli lab experiment =====
|
||||
|
||||
|
||||
!bnotice Typical pattern:
|
||||
The population grows faster and faster. "Why? Is there an underlying (general) mechanism":"http://www.zo.utexas.edu/courses/Thoc/PopGrowth.html"?
|
||||
!enotice
|
||||
!bblock
|
||||
o Cells divide after $T$ seconds on average (one generation)
|
||||
o $2N$ celles divide into twice as many new cells $\Delta N$ in a time
|
||||
interval $\Delta t$ as $N$ cells would: $\Delta N \propto N$
|
||||
o $N$ cells result in twice as many new individuals $\Delta N$ in
|
||||
time $2\Delta t$ as in time $\Delta t$: $\Delta N \propto\Delta t$
|
||||
o Same proportionality wrt death (repeat reasoning)
|
||||
o Proposed model: $\Delta N = b\Delta t N - d\Delta tN$ for some unknown
|
||||
constants $b$ (births) and $d$ (deaths)
|
||||
o Describe evolution in discrete time: $t_n=n\Delta t$
|
||||
o Program-friendly notation: $N$ at $t_n$ is $N^n$
|
||||
o Math model: $N^{n+1} = N^n + r\Delta t\, N$ (with $\ r=b-d$)
|
||||
o Program model: `N[n+1] = N[n] + r*dt*N[n]`
|
||||
!eblock
|
||||
|
||||
!split
|
||||
===== The program =====
|
||||
|
||||
!bblock
|
||||
Let us solve the difference equation in as simple way as possible,
|
||||
just to train some programming: $r=1.5$, $N^0=1$, $\Delta t=0.5$
|
||||
|
||||
@@@CODE ../../Programs/LifeScience/diffeq.py
|
||||
!eblock
|
||||
|
||||
% if FORMAT != 'ipynb':
|
||||
!split
|
||||
===== The output =====
|
||||
|
||||
!bc
|
||||
N[1]=1.2
|
||||
N[2]=1.6
|
||||
N[3]=2.0
|
||||
N[4]=2.4
|
||||
N[5]=3.1
|
||||
N[6]=3.8
|
||||
N[7]=4.8
|
||||
N[8]=6.0
|
||||
N[9]=7.5
|
||||
N[10]=9.3
|
||||
N[11]=11.6
|
||||
N[12]=14.6
|
||||
N[13]=18.2
|
||||
N[14]=22.7
|
||||
N[15]=28.4
|
||||
N[16]=35.5
|
||||
N[17]=44.4
|
||||
N[18]=55.5
|
||||
N[19]=69.4
|
||||
N[20]=86.7
|
||||
!ec
|
||||
% endif
|
||||
|
||||
!split
|
||||
===== Parameter estimation =====
|
||||
|
||||
!bblock
|
||||
* We do not know $r$
|
||||
* How can we estimate $r$ from data?
|
||||
|
||||
We can use the difference equation with the experimental data
|
||||
|
||||
!bt
|
||||
\[ N^{n+1} = N^n + r\Delta t N^n\]
|
||||
!et
|
||||
Say $N^{n+1}$ and $N^n$ are known from data, solve wrt $r$:
|
||||
|
||||
!bt
|
||||
\[ r = \frac{N^{n+1}-N^n}{N^n\Delta t} \]
|
||||
!et
|
||||
|
||||
Use experimental data in the fraction, say $t_1=600$, $t_2=1200$,
|
||||
$N^1=140$, $N^2=250$: $r=0.0013$.
|
||||
!eblock
|
||||
|
||||
|
||||
!split
|
||||
===== A program relevant for the biological problem =====
|
||||
|
||||
# exact r = 0.000694
|
||||
|
||||
!bblock
|
||||
@@@CODE ../../Programs/LifeScience/ecoli.py
|
||||
|
||||
Change `r` in the program and play around to make a better fit!
|
||||
!eblock
|
||||
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Simulating financial transcations =====
|
||||
|
||||
The aim here is to simulate financial transactions among financial agents
|
||||
using Monte Carlo methods. The final goal is to extract a distribution of income as function
|
||||
of the income $m$. From Pareto's work ("V.~Pareto, 1897":"http://www.institutcoppet.org/2012/05/08/cours-deconomie-politique-1896-de-vilfredo-pareto") it is known from empirical studies
|
||||
that the higher end of the distribution of money follows a distribution
|
||||
!bt
|
||||
\[
|
||||
w_m\propto m^{-1-\alpha},
|
||||
\]
|
||||
!et
|
||||
with $\alpha\in [1,2]$. We will here follow the analysis made by "Patriarca and collaborators":"http://www.sciencedirect.com/science/article/pii/S0378437104004327".
|
||||
|
||||
Here we will study numerically the relation between the micro-dynamic relations among financial
|
||||
agents and the resulting macroscopic money distribution.
|
||||
|
||||
We assume we have $N$ agents that exchange money in pairs $(i,j)$. We assume also that all agents
|
||||
start with the same amount of money $m_0 > 0$. At a given 'time step', we choose randomly a pair
|
||||
of agents $(i,j)$ and let a transaction take place. This means that agent $i$'s money $m_i$ changes
|
||||
to $m_i'$ and similarly we have $m_j\rightarrow m_j'$.
|
||||
Money is conserved during a transaction, meaning that
|
||||
!bt
|
||||
\begin{equation}
|
||||
m_i+m_j=m_i'+m_j'.
|
||||
label{eq:conserve}
|
||||
\end{equation}
|
||||
!et
|
||||
The change is done via a random reassignement (a random number) $\epsilon$, meaning that
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
m_i' = \epsilon(m_i+m_j),
|
||||
\end{equation*}
|
||||
!et
|
||||
leading to
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
m_j'= (1-\epsilon)(m_i+m_j).
|
||||
\end{equation*}
|
||||
!et
|
||||
The number $\epsilon$ is extracted from a uniform distribution.
|
||||
In this simple model, no agents are left with a debt, that is $m\ge 0$.
|
||||
Due to the conservation law above, one can show that the system relaxes toward an equilibrium
|
||||
state given by a Gibbs distribution
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
w_m=\beta \exp{(-\beta m)},
|
||||
\end{equation*}
|
||||
!et
|
||||
with
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
\beta = \frac{1}{\langle m\rangle},
|
||||
\end{equation*}
|
||||
!et
|
||||
and $\langle m\rangle=\sum_i m_i/N=m_0$, the average money.
|
||||
It means that after equilibrium has been reached that the majority of agents is left with a small
|
||||
number of money, while the number of richest agents, those with $m$ larger than a specific value $m'$,
|
||||
exponentially decreases with $m'$.
|
||||
|
||||
We assume that we have $N=500$ agents. In each simulation, we need a sufficiently large number of transactions, say $10^7$. Our aim is find the final equilibrium distribution $w_m$. In order to do that we would need
|
||||
several runs of the above simulations, at least $10^3-10^4$ runs (experiments).
|
||||
|
||||
|
||||
|
||||
|
||||
=== Simulation of Transactions ===
|
||||
Our task is to first set up an algorithm which simulates the above transactions with an initial
|
||||
amount $m_0$.
|
||||
The challenge here is to figure out a Monte Carlo simulation based on the
|
||||
above equations.
|
||||
You will in particular need to make an algorithm which sets up a histogram as function of $m$.
|
||||
This histogram contains the number of times a value $m$ is registered and represents
|
||||
$w_m\Delta m$. You will need to set up a value for the interval $\Delta m$ (typically $0.01-0.05$).
|
||||
That means you need to account for the number of times you register an income in the interval
|
||||
$m,m+\Delta m$. The number of times you register this income, represents the value that enters the histogram.
|
||||
You will also need to find a criterion for when the equilibrium situation has been reached.
|
||||
|
||||
!bc pycod
|
||||
#!/usr/bin/env python
|
||||
import numpy as np
|
||||
import matplotlib.mlab as mlab
|
||||
import matplotlib.pyplot as plt
|
||||
import random
|
||||
|
||||
# initialize the rng with a seed
|
||||
random.seed()
|
||||
# Hard coding of input parameters
|
||||
Agents = 500
|
||||
MCcounts = 1000
|
||||
Transactions = 100000
|
||||
startMoney = 1.0
|
||||
Lambda = 0.0
|
||||
FinancialAgents = startMoney*np.ones(Agents)
|
||||
for i in range (1, MCcounts, 1):
|
||||
for j in range (1, Transactions, 1):
|
||||
agent_i = int(Agents*random.random())
|
||||
agent_j = int(Agents*random.random())
|
||||
epsilon = random.random()
|
||||
if agent_i != agent_j:
|
||||
m1 = Lambda*FinancialAgents[agent_i] + (1-Lambda)*epsilon*(FinancialAgents[agent_i] + FinancialAgents[agent_j])
|
||||
m2 = Lambda*FinancialAgents[agent_j] + (1-Lambda)*(1-epsilon)*(FinancialAgents[agent_i] + FinancialAgents[agent_j])
|
||||
FinancialAgents[agent_i] = m1
|
||||
FinancialAgents[agent_j] = m2
|
||||
|
||||
# the histogram of the data
|
||||
n, bins, patches = plt.hist(FinancialAgents, 50, facecolor='green')
|
||||
|
||||
plt.xlabel('$x$')
|
||||
plt.ylabel('Distribution of wealth')
|
||||
plt.title(r'Money')
|
||||
plt.axis([0, 10, 0, 500])
|
||||
plt.grid(True)
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
We can then change our model to allow for a saving criterion, meaning that the agents save
|
||||
a fraction $\lambda$ of the money they have before the transaction is made. The final distribution will then no longer be given by Gibbs distribution. It could also include a taxation on financial transactions.
|
||||
|
||||
The conservation law of Eq. (ref{eq:conserve}) holds, but the money to be shared in a transaction between
|
||||
agent $i$ and agent $j$ is now $(1-\lambda)(m_i+m_j)$. This means that we have
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
m_i' = \lambda m_i+\epsilon(1-\lambda)(m_i+m_j),
|
||||
\end{equation*}
|
||||
!et
|
||||
and
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
m_j' = \lambda m_j+(1-\epsilon)(1-\lambda)(m_i+m_j),
|
||||
\end{equation*}
|
||||
!et
|
||||
which can be written as
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
m_i'=m_i+\delta m
|
||||
\end{equation*}
|
||||
!et
|
||||
and
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
m_j'=m_j-\delta m,
|
||||
\end{equation*}
|
||||
!et
|
||||
with
|
||||
|
||||
!bt
|
||||
\begin{equation*}
|
||||
\delta m=(1-\lambda)(\epsilon m_j-(1-\epsilon)m_i),
|
||||
\end{equation*}
|
||||
!et
|
||||
showing how money is conserved during a transaction.
|
||||
Select values of $\lambda =0.25,0.5$ and $\lambda=0.9$ and try to extract the corresponding
|
||||
equilibrium distributions and compare these with the Gibbs distribution. Comment your results.
|
||||
Extract a parametrization of the above curves, see for example "Patriarca and collaborators":"http://www.sciencedirect.com/science/article/pii/S0378437104004327" and see if you can parametrize the high-end tails of the distributions in terms of power laws. Comment your results.
|
||||
|
||||
In the studies above the agents were selected randomly, irrespective of whether we allowed for
|
||||
saving or not during a transaction. What is often observed is that various agents tend to make preferences for for whom to interact with. We will now study the evolution of the distribution of wealth $w_m$ by assuming that there is a likelihood
|
||||
!bt
|
||||
\[
|
||||
p_{ij} \propto \vert m_i-m_j\vert^{-\alpha},
|
||||
\]
|
||||
!et
|
||||
for an interaction between agents $i$ and $j$ with respective wealths $m_i$ and $m_j$. The parameter $\alpha > 0$. For $\alpha=0$ we recover our model from part 5a).
|
||||
Perform the same analysis as previously with $N=500$ as well as with $N=1000$ agents and study the distribution of wealth for $\alpha =0.5$, $\alpha =1.0$, $\alpha =1.5$ and $\alpha =2.0$.
|
||||
You should try to reproduce Figure 1 of "Goswami and Sen":"http://www.sciencedirect.com/science/article/pii/S0378437114006967".
|
||||
Extract the tail of the distribution and see if it follows a Pareto distribution
|
||||
!bt
|
||||
\[
|
||||
w_m\propto m^{-1-\alpha}.
|
||||
\]
|
||||
!et
|
||||
What happens if $\alpha \gg 1$?
|
||||
|
||||
Perform the analysis with and without a saving $\lambda$ on each transaction and comment your results.
|
||||
We add to the previous probability the possibility that two agents who interact have performed similar transactions earlier. That is, in addition to being financially close, we assume that the likelihood for interacting increases if two agents have interacted earlier.
|
||||
We add this feature by modifying the previous likelihood to
|
||||
!bt
|
||||
\[
|
||||
p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}\left(c_{ij}+1\right)^{\gamma},
|
||||
\]
|
||||
!et
|
||||
where $c_{ij}$ represents the number of previous interactions that have taken place between $i$ and $j$. The factor $1$ is added in order to ensure that if they have not interacted earlier they can still interact. Perform similar studies as above with $N=1000$, $\alpha=1.0$ and $\alpha=2.0$ using $\gamma = 0.0, 1.0, 2.0, 3.0$ and $4.0$. Plot the wealth distributions for these cases and try to extract eventual power law tails with and without a saving $\lambda$ in each transaction. Comment your results and compare them with figures 5 and 6 of "Goswami and Sen":"http://www.sciencedirect.com/science/article/pii/S0378437114006967".
|
||||
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 32 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 52 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 30 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 34 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 40 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 29 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 221 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 5.3 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 8.1 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 16 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 9.5 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 7.0 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 11 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 13 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 156 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 8.7 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 16 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 69 KiB |
@@ -44,7 +44,7 @@ system doconce split_html $html.html --method=space10
|
||||
# Bootstrap style
|
||||
html=${name}-bs
|
||||
system doconce format html $name --html_style=bootstrap --pygments_html_style=default --html_admon=bootstrap_panel --html_output=$html $opt
|
||||
system doconce split_html $html.html --method=split --pagination --nav_button=bottom
|
||||
#system doconce split_html $html.html --method=split --pagination --nav_button=bottom
|
||||
|
||||
# IPython notebook
|
||||
system doconce format ipynb $name $opt
|
||||
|
||||
@@ -5,15 +5,19 @@ DATE: today
|
||||
!split
|
||||
===== What is Machine Learning? =====
|
||||
|
||||
Machine learning is the science of giving computers the ability to learn without being explicitly programmed.
|
||||
The idea is that there exist generic algorithms which can be used to find patterns in a broad class of data sets without
|
||||
having to write code specifically for each problem. The algorithm will build its own logic based on the data.
|
||||
Machine learning is the science of giving computers the ability to
|
||||
learn without being explicitly programmed. The idea is that there
|
||||
exist generic algorithms which can be used to find patterns in a broad
|
||||
class of data sets without having to write code specifically for each
|
||||
problem. The algorithm will build its own logic based on the data.
|
||||
|
||||
Machine learning is a subfield of computer science, and is closely related to computational statistics.
|
||||
It evolved from the study of pattern recognition in artificial intelligence (AI) research, and has made contributions to
|
||||
AI tasks like computer vision, natural language processing
|
||||
and speech recognition. It has also, especially in later years,
|
||||
found applications in a wide variety of other areas, including bioinformatics, economy, physics, finance and marketing.
|
||||
Machine learning is a subfield of computer science, and is closely
|
||||
related to computational statistics. It evolved from the study of
|
||||
pattern recognition in artificial intelligence (AI) research, and has
|
||||
made contributions to AI tasks like computer vision, natural language
|
||||
processing and speech recognition. It has also, especially in later
|
||||
years, found applications in a wide variety of other areas, including
|
||||
bioinformatics, economy, physics, finance and marketing.
|
||||
|
||||
!split
|
||||
===== Types of Machine Learning =====
|
||||
@@ -39,14 +43,18 @@ Some of the most common tasks are:
|
||||
|
||||
!split
|
||||
===== Artificial neurons =====
|
||||
The field of artificial neural networks has a long history of development, and is closely connected with
|
||||
the advancement of computer science and computers in general. A model of artificial neurons
|
||||
was first developed by McCulloch and Pitts in 1943 to study signal processing in the brain and
|
||||
has later been refined by others. The general idea is to mimic neural networks in the human brain, which
|
||||
is composed of billions of neurons that communicate with each other by sending electrical signals.
|
||||
Each neuron accumulates its incoming signals,
|
||||
which must exceed an activation threshold to yield an output. If the threshold is not overcome, the neuron
|
||||
remains inactive, i.e. has zero output.
|
||||
|
||||
The field of artificial neural networks has a long history of
|
||||
development, and is closely connected with the advancement of computer
|
||||
science and computers in general. A model of artificial neurons was
|
||||
first developed by McCulloch and Pitts in 1943 to study signal
|
||||
processing in the brain and has later been refined by others. The
|
||||
general idea is to mimic neural networks in the human brain, which is
|
||||
composed of billions of neurons that communicate with each other by
|
||||
sending electrical signals. Each neuron accumulates its incoming
|
||||
signals, which must exceed an activation threshold to yield an
|
||||
output. If the threshold is not overcome, the neuron remains inactive,
|
||||
i.e. has zero output.
|
||||
|
||||
This behaviour has inspired a simple mathematical model for an artificial neuron.
|
||||
|
||||
@@ -59,21 +67,57 @@ This behaviour has inspired a simple mathematical model for an artificial neuron
|
||||
Here, the output $y$ of the neuron is the value of its activation function, which have as input
|
||||
a weighted sum of signals $x_i, \dots ,x_n$ received by $n$ other neurons.
|
||||
|
||||
Conceptually, it is helpful to divide neural networks into four
|
||||
categories:
|
||||
o general purpose neural networks for supervised learning,
|
||||
o neural networks designed specifically for image processing, the most prominent example of this class being Convolutional Neural Networks (CNNs),
|
||||
o neural networks for sequential data such as Recurrent Neural Networks (RNNs), and
|
||||
o neural networks for unsupervised learning such as Deep Boltzmann Machines.
|
||||
|
||||
|
||||
In physics, DNNs and CNNs have already found numerous applications. In
|
||||
statistical physics, they have been applied to detect phase
|
||||
transitions in 2D Ising and Potts models, lattice gauge theories, and
|
||||
different phases of polymers.
|
||||
Deep learning has also found interesting applications in quantum
|
||||
physics. Various quantum phase transitions can be detected and studied
|
||||
using DNNs and CNNs, including the transverse-field Ising model,
|
||||
topological phases, and even non-equilibrium many-body
|
||||
localization. Representing quantum states as DNNs quantum state
|
||||
tomography are among some of the impressive
|
||||
achievements to reveal the potential of DNNs to facilitate the study
|
||||
of quantum systems.
|
||||
|
||||
In quantum information theory, it has been shown that one can perform
|
||||
gate decompositions with the help of neural. In lattice quantum chromodynamics,
|
||||
DNNs have been used to learn action parameters in regions of parameter
|
||||
space where PCA fails. Last but not least,
|
||||
DNNs also found place in the study of quantum, and in scattering theory to learn
|
||||
$s$-wave scattering length of potentials.
|
||||
|
||||
!split
|
||||
===== Neural network types =====
|
||||
|
||||
An artificial neural network (NN), is a computational model that consists of layers of connected neurons, or *nodes*.
|
||||
It is supposed to mimic a biological nervous system by letting each neuron interact with other neurons
|
||||
by sending signals in the form of mathematical functions between layers.
|
||||
A wide variety of different NNs have
|
||||
been developed, but most of them consist of an input layer, an output layer and eventual layers in-between, called
|
||||
*hidden layers*. All layers can contain an arbitrary number of nodes, and each connection between two nodes
|
||||
is associated with a weight variable.
|
||||
An artificial neural network (NN), is a computational model that
|
||||
consists of layers of connected neurons, or *nodes*. It is supposed
|
||||
to mimic a biological nervous system by letting each neuron interact
|
||||
with other neurons by sending signals in the form of mathematical
|
||||
functions between layers. A wide variety of different NNs have been
|
||||
developed, but most of them consist of an input layer, an output layer
|
||||
and eventual layers in-between, called *hidden layers*. All layers can
|
||||
contain an arbitrary number of nodes, and each connection between two
|
||||
nodes is associated with a weight variable.
|
||||
|
||||
Neural networks (also called neural nets) are neural-inspired
|
||||
nonlinear models for supervised learning. As we will see, neural nets
|
||||
can be viewed as natural, more powerful extensions of supervised
|
||||
learning methods such as linear and logistic regression and soft-max
|
||||
methods.
|
||||
|
||||
|
||||
!split
|
||||
===== Feed-forward neural networks =====
|
||||
|
||||
The feed-forward neural network (FFNN) was the first and simplest type of NN devised. In this network,
|
||||
the information moves in only one direction: forward through the layers.
|
||||
|
||||
@@ -351,3 +395,5 @@ performs better than the sigmoid for training MLPs.
|
||||
has become the most popular for *deep neural networks*
|
||||
|
||||
|
||||
!split
|
||||
===
|
||||
|
||||
@@ -1,16 +1,14 @@
|
||||
TITLE: Project on Machine Learning
|
||||
AUTHOR: "Data Analysis and Machine Learning FYS-MAT3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html" {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo, Norway
|
||||
DATE: Fall semester 2017
|
||||
AUTHOR: "Data Analysis and Machine Learning FYS-STK3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html" {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo, Norway
|
||||
DATE: May 2018
|
||||
|
||||
|
||||
===== Using results from Monte Carlo models for machine learning =====
|
||||
===== Machine learning (ML) approaches to data from Ising model calculations =====
|
||||
|
||||
=== Introduction ===
|
||||
|
||||
The aim of this project is to use an already developed Monte Carlo program (either the ising Model or a variational Monte Carlo code) to produce, in case of the "Ising model":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/IsingModel", the energy as function of temperature. For the "variational Monte Carlo calculation":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/VMC2Electrons/" of interacting electrons in an oscilaltor trap, the data are represented by the ground state energies as functions of the variational parameters.
|
||||
|
||||
In its simplest form
|
||||
the energy of the Ising model is expressed as, without an externally applied magnetic field,
|
||||
The aim of this project is to use an already developed Monte Carlo program for the one-dimensional and two-dimensional "Ising model":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/IsingModel", in order to produce the spin configurations for a series of energies $E_i$ (10000 in total) for a system of $L=40$ spins in one dimension and $L=40\times 40$ in two dimensions at three different temperatures.
|
||||
In its simplest form the energy of the Ising model is expressed as, without an externally applied magnetic field,
|
||||
!bt
|
||||
\[
|
||||
E=-J\sum_{< kl >}^{N}s_ks_l
|
||||
@@ -22,222 +20,86 @@ constant expressing the strength of the interaction between
|
||||
neighboring spins. The symbol $<kl>$ indicates that we sum over
|
||||
nearest neighbors only. We will assume that we have a ferromagnetic
|
||||
ordering, viz $J> 0$. We will use periodic boundary conditions and
|
||||
the Metropolis algorithm only. Alternatively we can use the supplied variational Monte Carlo program which solves Schroedinger's equation for two interacting electrons in a harmonic oscillator trap. Both codes can be found at the webpage of the "course under programs":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/".
|
||||
the Metropolis algorithm only. The spins take values $-1$ and $+1$ only.
|
||||
|
||||
|
||||
We will use the Ising model to generate our training data and will focus mainly on supervised training. We will follow closely the recent article of "Mehta et al, arXiv 1803.08823":"https://arxiv.org/abs/1803.08823". This article stands out as an excellent review on machine learning (ML) algorithms applied to typical physics problems. The added benefit is that each figure and model presented in "this article is accompanied by its jupyter notebook":"https://physics.bu.edu/~pankajm/MLnotebooks.html". This means that we can start using these and compare with our own results. In case you wish to use their data for the Ising model, their data can be downloaded from the same link which lists to the jupyter notebooks. See also at the end of the project description for more information on how to install various Python packages.
|
||||
|
||||
|
||||
|
||||
With the abovementioned configurations we will determine, using first various
|
||||
regression methods, the value of the coupling constant for the energy
|
||||
of the one-dimensional Ising model. Thereafter, we will use the
|
||||
two-dimensional data, but now computed at different temperatures, in
|
||||
order to classify the phase of the Ising model. Below the critical
|
||||
temperature, the system will be in a so-called ferromagnetic
|
||||
phase. Close to the critical temperature, the final magnetization becomes smaller and smaller in absolute value
|
||||
while above the critical temperature,
|
||||
the net magnetization is zero. This classification case, that is the
|
||||
two-dimensional Ising model, will be studied using logistic regression, a _random forest_
|
||||
algorithm and deep neural networks.
|
||||
|
||||
You should try to program at least one of these methods yourself (choose the one you prefer).
|
||||
Feel free to use the notebooks to benchmark your code. If you wish to write your own C++ or Fortran program for say a simple neural network model, please feel free to do so.
|
||||
You can then benchmark your results against the above jupyter notebooks. More information can also be found at the link for the lecture notes of "FYS-STK4155":"https://compphysics.github.io/MachineLearning/doc/web/course.html".
|
||||
|
||||
|
||||
We recommend that you form groups of 2-3 students and try to
|
||||
collaborate on the notebooks, develop your own software and discuss
|
||||
the final presentations. You can collaborate on all these topics. The
|
||||
final presentation should include an overview of popular machine
|
||||
learning algorithms as introduction and motivation. Thereafter you
|
||||
discuss the explicit Ising model data and how you have implemented the
|
||||
ML algorithms discussed here, discuss their pros and cons and try to
|
||||
develop your own code for at least one of these algorithms. You are
|
||||
encouraged to use the abovementioned notebooks as starting point and
|
||||
guidance. The duration of your presentation should at most be 30
|
||||
mins. Allow for approximately 15 mins for discussions and questions.
|
||||
|
||||
|
||||
|
||||
=== Part a): Producing the data ===
|
||||
|
||||
If we opt for the Ising model code, we need to generate for every span over the lattice and output of the total energy and the magnetic moment (if we wish to study this quantity as well), that is compute and write to file as function of temperature
|
||||
$\langle E\rangle$ and $\langle \vert M\vert \rangle$.
|
||||
We will use a fixed lattice siste of $L\times L = 40 \times 40$. Make sure the calculations have been equilibrated and compute these expectation values for
|
||||
for $T\in [1.0,3.0]$ with a step in
|
||||
temperature $\Delta T=0.1$ or smaller. This are the data you will use to estimate the standard deviation in the next part of the project.
|
||||
You should keep a file for the different temperatures as you will need to compute the standard deviation for each temperature.
|
||||
You can use the Ising model data from the article of Mehta *et al.*, or generate your own data.
|
||||
If you opt for using your own Ising model code, you need to generate $10000$ energy configurations with their spin orientations after the system has reached its most likely state. These energies and their corresponding spin orientations
|
||||
represent then your data.
|
||||
We will use a fixed lattice of $L\times L = 40 \times 40$ spins in two dimensions and $L=40$ spins in one dimension.
|
||||
Make sure the calculations have been equilibrated. For the two-dimensional system, compute the configurations
|
||||
for three values of the temperature, namely $T=0.75$ (ordered phase), $T=2.3$ (near the critical point) and $T=4.0$ (disordered phase).
|
||||
For the one-dimensional system it suffices to compute the various configurations for one temperature only, say $T=2.0$.
|
||||
These are the data you will use to study different ML algorithms.
|
||||
We generate our data with $J=1$.
|
||||
|
||||
=== Part b): Estimating the coupling constant of the one-dimensional Ising model ===
|
||||
|
||||
We start with the one-dimensional Ising model and use the data we have generated with $J=1$. Use linear regression, Lasso and Ridge regression as described section 6 and in Notebook 4 of "Mehta *et al.*":"https://physics.bu.edu/~pankajm/ML-Notebooks/HTML/NB_CVI-linreg_ising.html". Discuss the methods and how they perform in computing the coupling constant $J$. Give a critical analysis and discuss how to evaluate the *cost function*. You should feel free to write your own code, see also
|
||||
the lecture notes of "FYS-STK4155":"https://compphysics.github.io/MachineLearning/doc/web/course.html", in particular te material on least square methods. You can use scikit-learn to perform these analyses. See below for instruction on how to install scikit-learn.
|
||||
|
||||
=== Part c): Determine the phase of the two-dimensional Ising model ===
|
||||
|
||||
We switch now to binary classification methods and use logistic regression to define the phases of the Ising model.
|
||||
Use described section 7 and in Notebook 6 of "Mehta *et al.*":"https://physics.bu.edu/~pankajm/ML-Notebooks/HTML/NB_CVII-logreg_ising.html". Discuss the methods and how they perform. Give a critical analysis and discuss how to evaluate the *cost function*. You should feel free to write your own code. Use thereafter the *random forests* algorithm to classify the same phases as done with logistic regression and discuss the pros and cons of these methods. For *random forests* you can use "notebook 9":"https://physics.bu.edu/~pankajm/ML-Notebooks/HTML/NB_CVIII-randomforests_ising.html" of Mehta *et al.*
|
||||
|
||||
You can use scikit-learn to perform these analyses. See below for instruction on how to install scikit-learn.
|
||||
|
||||
=== Part d): Classifying the Ising model phase using neural networks ===
|
||||
|
||||
We end the classification problem of the phases of the Ising model by employing the algorithm for so-called feed-forward deep neural networks (see section 9 of Mehta *et al.*). The method is described in "notebook 12":"https://physics.bu.edu/~pankajm/ML-Notebooks/HTML/NB_CIX-DNN_ising_TFlow.html".
|
||||
|
||||
You can use tensorflow to perform these analyses. See below for instruction on how to install tensorflow.
|
||||
|
||||
|
||||
Alternatively you can run the variational Monte Carlo program for two interacting electrons confined to move in a harmonic oscillator trap. Here you will produce the variational expectaction value of the energy as a function of the variational parameters $\alpha$ and $\beta$. You could single out only a set of such parameters after you have found the minimum energy.
|
||||
=== Part e): Summary ===
|
||||
|
||||
These two cases form then our training data which we will use for our estimates of the standard deviation.
|
||||
|
||||
|
||||
=== Part b): Estimating the standard deviation ===
|
||||
|
||||
This part deals with widely used resampling methods to find the standard deviation on our data sets. The methods we can use are the Blocking method, the Bootstrap method and the Jackknife method. The latter two are examples of so-called resampling methods. "The functions provided under the program folder":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/Sampling" encode all three methods. The typical situation is to generate data using a C++ or Fortan code, as in part a) above. These data are then used by a post-analysis program in order to perform a statistical analysis.
|
||||
|
||||
In case of the Ising model, the final product will be a table of the expectation values of the energy (or the magnetic moment) as function of temperature and with a proper standard deviation. These data will then enter our analysis in parts c-e). We will call these data the training data.
|
||||
|
||||
If you opt for the variational Monte Carlo program, you need to make a table of the ground state energy as function of various variational parameters (you can limit yourself to $\alpha$ only) with proper standard deviation estimations.
|
||||
|
||||
Your task here is thus to use the data from part a) and generate the standard deviation using one (or more) of the methods discussed above. In general, the Bootstrap method is the most widely used one.
|
||||
|
||||
These data will then be used in our next step, where we will try to obtain a functional approximation to our data sets.
|
||||
|
||||
=== Part c): Fitting the data using regression analysis and other methods ===
|
||||
|
||||
With the data from part b), with a proper determination of the standard deviation, the task here is to find a fit to the data using regression methods. Here we will focus on the least square method with a straight line fit first and thereafter using singular value decomposition to fit a function that reproduces the data.
|
||||
|
||||
Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of $y$ which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree $n-1$ with $n$ points, that is
|
||||
!bt
|
||||
\[
|
||||
y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_i x_i^j+\epsilon_i,
|
||||
\]
|
||||
!et
|
||||
where $\epsilon_i$ is the error in our approximation.
|
||||
|
||||
In order to find the optimal parameters $\beta_i$ instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values $y_i$ (which represent hopefully the exact values) and the parametrized values $\tilde{y}_i$, namely
|
||||
!bt
|
||||
\[
|
||||
Q(\hat{\beta})=\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\left(\hat{y}-\hat{\tilde{y}}\right)^T\left(\hat{y}-\hat{\tilde{y}}\right),
|
||||
\]
|
||||
!et
|
||||
or using the matrix $\hat{X}$ as
|
||||
!bt
|
||||
\[
|
||||
Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right).
|
||||
\]
|
||||
!et
|
||||
|
||||
The function
|
||||
!bt
|
||||
\[
|
||||
Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right),
|
||||
\]
|
||||
!et
|
||||
can be linked to the variance of the quantity $y_i$ if we interpret the latter as the mean value of for example a numerical experiment. When linking below with the maximum likelihood approach below, we will indeed interpret $y_i$ as a mean value
|
||||
!bt
|
||||
\[
|
||||
y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
|
||||
\]
|
||||
!et
|
||||
where $\langle y_i \rangle$ is the mean value. Keep in mind also that till now we have treated $y_i$ as the exact value. Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.
|
||||
|
||||
Introducing the standard deviation $\sigma_i$ for each measurement $y_i$, we define now the $\chi^2$ function as
|
||||
!bt
|
||||
\[
|
||||
\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right),
|
||||
\]
|
||||
!et
|
||||
where the matrix $\hat{\Sigma}$ is a diagonal matrix with $\sigma_i$ as matrix elements.
|
||||
|
||||
|
||||
|
||||
|
||||
In order to find the parameters $\beta_i$ we will then minimize the spread of $\chi^2(\hat{\beta})$ by requiring
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
|
||||
\]
|
||||
!et
|
||||
which results in
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
|
||||
\]
|
||||
!et
|
||||
or in a matrix-vector form as
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right).
|
||||
\]
|
||||
!et
|
||||
where we have defined the matrix $\hat{A} =\hat{X}/\hat{\Sigma}$ with matrix elements $a_{ij} = x_{ij}/\sigma_i$ and the vector $\hat{b}$ with elements $b_i = y_i/\sigma_i$.
|
||||
|
||||
|
||||
We can rewrite
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right),
|
||||
\]
|
||||
!et
|
||||
as
|
||||
!bt
|
||||
\[
|
||||
\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta},
|
||||
\]
|
||||
!et
|
||||
and if the matrix $\hat{A}^T\hat{A}$ is invertible we have the solution
|
||||
!bt
|
||||
\[
|
||||
\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}.
|
||||
\]
|
||||
!et
|
||||
If we then introduce the matrix
|
||||
!bt
|
||||
\[
|
||||
\hat{H} = \hat{A}^T\hat{A},
|
||||
\]
|
||||
!et
|
||||
we have then the following expression for the parameters $\beta_j$ (the matrix elements of $\hat{H}$ are $h_{ij}$)
|
||||
!bt
|
||||
\[
|
||||
\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik}
|
||||
\]
|
||||
!et
|
||||
We state without proof the expression for the uncertainty in the parameters $\beta_j$ as
|
||||
!bt
|
||||
\[
|
||||
\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
|
||||
\]
|
||||
!et
|
||||
resulting in
|
||||
!bt
|
||||
\[
|
||||
\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}!
|
||||
\]
|
||||
!et
|
||||
|
||||
|
||||
The first step here is to approximate the function $y$ with a first-order polynomial, that is we write
|
||||
!bt
|
||||
\[
|
||||
y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
|
||||
\]
|
||||
!et
|
||||
By computing the derivatives of $\chi^2$ with respect to $\beta_0$ and $\beta_1$ show that these are given by
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
|
||||
\]
|
||||
!et
|
||||
and
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
|
||||
\]
|
||||
!et
|
||||
Define then
|
||||
!bt
|
||||
\[
|
||||
\gamma = \sum_{i=0}^{1}\frac{1}{\sigma_i^2},
|
||||
\]
|
||||
!et
|
||||
|
||||
!bt
|
||||
\[
|
||||
\gamma_x = \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2},
|
||||
\]
|
||||
!et
|
||||
!bt
|
||||
\[
|
||||
\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right),
|
||||
\]
|
||||
!et
|
||||
!bt
|
||||
\[
|
||||
\gamma_{xx} = \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2},
|
||||
\]
|
||||
!et
|
||||
!bt
|
||||
\[
|
||||
\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2},
|
||||
\]
|
||||
!et
|
||||
and show that
|
||||
!bt
|
||||
\[
|
||||
\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
|
||||
\]
|
||||
!et
|
||||
!bt
|
||||
\[
|
||||
\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
|
||||
\]
|
||||
!et
|
||||
Find these parameters for your data set and compare the fitted function $y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i$ and find the error in the parameters $\beta_0$ and $\beta_1$ as well.
|
||||
How does your parametrization compare with the input data $y_i$? Give a critical discussion of the method.
|
||||
|
||||
Develop now a program where you extract the parameters $\beta_i$ and their corresponding errors using the general Least Square method (LSM), with polynomials of degree $p=2$ and $p=3$. Which polyomials gives the best fit to your data? Here you should prepare an input file with the number of data points in your training data set, the corresponding standard deviation $\sigma_i$ as well as the degree $p$ of the polynomial you want to fit.
|
||||
|
||||
The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients $\beta_i$. A better approach is to use the Singular Value Decomposition (SVD) method discussed in the lecture notes on regression analysis. Here you should use the "SVD programs":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/SVD" and use the same input file you used for the general LSM. In the main functions you should specifiy which function you wish to approximate the date with. The example programs discussed above show you examples on how to use these codes as well as how to write the functional form of the polynomial or fitting function.
|
||||
|
||||
|
||||
=== Part d): Introducing Bayesian statistics ===
|
||||
More text to come
|
||||
|
||||
=== Part e): Studying the Ising model or the VMC results with Neural networks ===
|
||||
More text to come
|
||||
You should make a summary of the various methods and their pros and cons. For the final presentation, you should have at least coded one of these methods yourself and discussed the evaluation of the cost function.
|
||||
|
||||
===== Background literature =====
|
||||
|
||||
On Machine Learning we recommend strongly the article of Mehta *et al.*
|
||||
Textbooks on Machine Learning can be found at the "Github address of FYS-STK4155":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Textbooks", see in particular Marsland's text.
|
||||
|
||||
* "Mehta et al, arXiv 1803.08823":"https://arxiv.org/abs/1803.08823", *A high-bias, low-variance introduction to Machine Learning for physicists*, ArXiv:1803.08823.
|
||||
|
||||
If you wish to read more about the Ising model and statistical physics here are three suggestions.
|
||||
|
||||
* "M. Plischke and B. Bergersen":"http://www.worldscientific.com/worldscibooks/10.1142/5660", *Equilibrium Statistical Physics*, World Scientific, see chapters 5 and 6.
|
||||
@@ -271,10 +133,35 @@ project.
|
||||
* Try to establish a practice where you log your work at the computerlab. You may find such a logbook very handy at later stages in your work, especially when you don't properly remember what a previous test version of your program did. Here you could also record the time spent on solving the exercise, various algorithms you may have tested or other topics which you feel worthy of mentioning.
|
||||
|
||||
|
||||
===== Software and needed installations =====
|
||||
|
||||
If you have Python installed (we recommend Python3) and you feel pretty familiar with installing different packages,
|
||||
we recommend that you install the following Python packages via _pip_ as
|
||||
o pip install numpy scipy matplotlib ipython scikit-learn tensorflow sympy pandas pillow
|
||||
For Python3, replace _pip_ with _pip3_.
|
||||
|
||||
See below for a discussion of _tensorflow_ and _scikit-learn_.
|
||||
|
||||
For OSX users we recommend also, after having installed Xcode, to install _brew_. Brew allows
|
||||
for a seamless installation of additional software via for example
|
||||
o brew install python3
|
||||
|
||||
For Linux users, with its variety of distributions like for example the widely popular Ubuntu distribution
|
||||
you can use _pip_ as well and simply install Python as
|
||||
o sudo apt-get install python3 (or python for python2.7)
|
||||
etc etc.
|
||||
|
||||
If you don't want to install various Python packages with their dependencies separately, we recommend two widely used distrubutions which set up all relevant dependencies for Python, namely
|
||||
o "Anaconda":"https://docs.anaconda.com/" Anaconda is an open source distribution of the Python and R programming languages for large-scale data processing, predictive analytics, and scientific computing, that aims to simplify package management and deployment. Package versions are managed by the package management system _conda_
|
||||
o "Enthought canopy":"https://www.enthought.com/product/canopy/" is a Python distribution for scientific and analytic computing distribution and analysis environment, available for free and under a commercial license.
|
||||
|
||||
Popular software packages written in Python for ML are
|
||||
|
||||
* "Scikit-learn":"http://scikit-learn.org/stable/",
|
||||
* "Tensorflow":"https://www.tensorflow.org/",
|
||||
* "PyTorch":"http://pytorch.org/" and
|
||||
* "Keras":"https://keras.io/".
|
||||
These are all freely available at their respective GitHub sites. They
|
||||
encompass communities of developers in the thousands or more. And the number
|
||||
of code developers and contributors keeps increasing.
|
||||
|
||||
|
||||
@@ -0,0 +1,280 @@
|
||||
TITLE: Project on Machine Learning
|
||||
AUTHOR: "Data Analysis and Machine Learning FYS-MAT3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html" {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo, Norway
|
||||
DATE: Fall semester 2017
|
||||
|
||||
|
||||
===== Using results from Monte Carlo models for machine learning =====
|
||||
|
||||
=== Introduction ===
|
||||
|
||||
The aim of this project is to use an already developed Monte Carlo program (either the ising Model or a variational Monte Carlo code) to produce, in case of the "Ising model":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/IsingModel", the energy as function of temperature. For the "variational Monte Carlo calculation":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/VMC2Electrons/" of interacting electrons in an oscilaltor trap, the data are represented by the ground state energies as functions of the variational parameters.
|
||||
|
||||
In its simplest form
|
||||
the energy of the Ising model is expressed as, without an externally applied magnetic field,
|
||||
!bt
|
||||
\[
|
||||
E=-J\sum_{< kl >}^{N}s_ks_l
|
||||
\]
|
||||
!et
|
||||
with
|
||||
$s_k=\pm 1$. The quantity $N$ represents the total number of spins and $J$ is a coupling
|
||||
constant expressing the strength of the interaction between
|
||||
neighboring spins. The symbol $<kl>$ indicates that we sum over
|
||||
nearest neighbors only. We will assume that we have a ferromagnetic
|
||||
ordering, viz $J> 0$. We will use periodic boundary conditions and
|
||||
the Metropolis algorithm only. Alternatively we can use the supplied variational Monte Carlo program which solves Schroedinger's equation for two interacting electrons in a harmonic oscillator trap. Both codes can be found at the webpage of the "course under programs":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/".
|
||||
|
||||
=== Part a): Producing the data ===
|
||||
|
||||
If we opt for the Ising model code, we need to generate for every span over the lattice and output of the total energy and the magnetic moment (if we wish to study this quantity as well), that is compute and write to file as function of temperature
|
||||
$\langle E\rangle$ and $\langle \vert M\vert \rangle$.
|
||||
We will use a fixed lattice siste of $L\times L = 40 \times 40$. Make sure the calculations have been equilibrated and compute these expectation values for
|
||||
for $T\in [1.0,3.0]$ with a step in
|
||||
temperature $\Delta T=0.1$ or smaller. This are the data you will use to estimate the standard deviation in the next part of the project.
|
||||
You should keep a file for the different temperatures as you will need to compute the standard deviation for each temperature.
|
||||
|
||||
|
||||
Alternatively you can run the variational Monte Carlo program for two interacting electrons confined to move in a harmonic oscillator trap. Here you will produce the variational expectaction value of the energy as a function of the variational parameters $\alpha$ and $\beta$. You could single out only a set of such parameters after you have found the minimum energy.
|
||||
|
||||
These two cases form then our training data which we will use for our estimates of the standard deviation.
|
||||
|
||||
|
||||
=== Part b): Estimating the standard deviation ===
|
||||
|
||||
This part deals with widely used resampling methods to find the standard deviation on our data sets. The methods we can use are the Blocking method, the Bootstrap method and the Jackknife method. The latter two are examples of so-called resampling methods. "The functions provided under the program folder":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/Sampling" encode all three methods. The typical situation is to generate data using a C++ or Fortan code, as in part a) above. These data are then used by a post-analysis program in order to perform a statistical analysis.
|
||||
|
||||
In case of the Ising model, the final product will be a table of the expectation values of the energy (or the magnetic moment) as function of temperature and with a proper standard deviation. These data will then enter our analysis in parts c-e). We will call these data the training data.
|
||||
|
||||
If you opt for the variational Monte Carlo program, you need to make a table of the ground state energy as function of various variational parameters (you can limit yourself to $\alpha$ only) with proper standard deviation estimations.
|
||||
|
||||
Your task here is thus to use the data from part a) and generate the standard deviation using one (or more) of the methods discussed above. In general, the Bootstrap method is the most widely used one.
|
||||
|
||||
These data will then be used in our next step, where we will try to obtain a functional approximation to our data sets.
|
||||
|
||||
=== Part c): Fitting the data using regression analysis and other methods ===
|
||||
|
||||
With the data from part b), with a proper determination of the standard deviation, the task here is to find a fit to the data using regression methods. Here we will focus on the least square method with a straight line fit first and thereafter using singular value decomposition to fit a function that reproduces the data.
|
||||
|
||||
Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of $y$ which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree $n-1$ with $n$ points, that is
|
||||
!bt
|
||||
\[
|
||||
y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_i x_i^j+\epsilon_i,
|
||||
\]
|
||||
!et
|
||||
where $\epsilon_i$ is the error in our approximation.
|
||||
|
||||
In order to find the optimal parameters $\beta_i$ instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values $y_i$ (which represent hopefully the exact values) and the parametrized values $\tilde{y}_i$, namely
|
||||
!bt
|
||||
\[
|
||||
Q(\hat{\beta})=\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\left(\hat{y}-\hat{\tilde{y}}\right)^T\left(\hat{y}-\hat{\tilde{y}}\right),
|
||||
\]
|
||||
!et
|
||||
or using the matrix $\hat{X}$ as
|
||||
!bt
|
||||
\[
|
||||
Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right).
|
||||
\]
|
||||
!et
|
||||
|
||||
The function
|
||||
!bt
|
||||
\[
|
||||
Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right),
|
||||
\]
|
||||
!et
|
||||
can be linked to the variance of the quantity $y_i$ if we interpret the latter as the mean value of for example a numerical experiment. When linking below with the maximum likelihood approach below, we will indeed interpret $y_i$ as a mean value
|
||||
!bt
|
||||
\[
|
||||
y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
|
||||
\]
|
||||
!et
|
||||
where $\langle y_i \rangle$ is the mean value. Keep in mind also that till now we have treated $y_i$ as the exact value. Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.
|
||||
|
||||
Introducing the standard deviation $\sigma_i$ for each measurement $y_i$, we define now the $\chi^2$ function as
|
||||
!bt
|
||||
\[
|
||||
\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right),
|
||||
\]
|
||||
!et
|
||||
where the matrix $\hat{\Sigma}$ is a diagonal matrix with $\sigma_i$ as matrix elements.
|
||||
|
||||
|
||||
|
||||
|
||||
In order to find the parameters $\beta_i$ we will then minimize the spread of $\chi^2(\hat{\beta})$ by requiring
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
|
||||
\]
|
||||
!et
|
||||
which results in
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
|
||||
\]
|
||||
!et
|
||||
or in a matrix-vector form as
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right).
|
||||
\]
|
||||
!et
|
||||
where we have defined the matrix $\hat{A} =\hat{X}/\hat{\Sigma}$ with matrix elements $a_{ij} = x_{ij}/\sigma_i$ and the vector $\hat{b}$ with elements $b_i = y_i/\sigma_i$.
|
||||
|
||||
|
||||
We can rewrite
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right),
|
||||
\]
|
||||
!et
|
||||
as
|
||||
!bt
|
||||
\[
|
||||
\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta},
|
||||
\]
|
||||
!et
|
||||
and if the matrix $\hat{A}^T\hat{A}$ is invertible we have the solution
|
||||
!bt
|
||||
\[
|
||||
\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}.
|
||||
\]
|
||||
!et
|
||||
If we then introduce the matrix
|
||||
!bt
|
||||
\[
|
||||
\hat{H} = \hat{A}^T\hat{A},
|
||||
\]
|
||||
!et
|
||||
we have then the following expression for the parameters $\beta_j$ (the matrix elements of $\hat{H}$ are $h_{ij}$)
|
||||
!bt
|
||||
\[
|
||||
\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik}
|
||||
\]
|
||||
!et
|
||||
We state without proof the expression for the uncertainty in the parameters $\beta_j$ as
|
||||
!bt
|
||||
\[
|
||||
\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
|
||||
\]
|
||||
!et
|
||||
resulting in
|
||||
!bt
|
||||
\[
|
||||
\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}!
|
||||
\]
|
||||
!et
|
||||
|
||||
|
||||
The first step here is to approximate the function $y$ with a first-order polynomial, that is we write
|
||||
!bt
|
||||
\[
|
||||
y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
|
||||
\]
|
||||
!et
|
||||
By computing the derivatives of $\chi^2$ with respect to $\beta_0$ and $\beta_1$ show that these are given by
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
|
||||
\]
|
||||
!et
|
||||
and
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
|
||||
\]
|
||||
!et
|
||||
Define then
|
||||
!bt
|
||||
\[
|
||||
\gamma = \sum_{i=0}^{1}\frac{1}{\sigma_i^2},
|
||||
\]
|
||||
!et
|
||||
|
||||
!bt
|
||||
\[
|
||||
\gamma_x = \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2},
|
||||
\]
|
||||
!et
|
||||
!bt
|
||||
\[
|
||||
\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right),
|
||||
\]
|
||||
!et
|
||||
!bt
|
||||
\[
|
||||
\gamma_{xx} = \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2},
|
||||
\]
|
||||
!et
|
||||
!bt
|
||||
\[
|
||||
\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2},
|
||||
\]
|
||||
!et
|
||||
and show that
|
||||
!bt
|
||||
\[
|
||||
\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
|
||||
\]
|
||||
!et
|
||||
!bt
|
||||
\[
|
||||
\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
|
||||
\]
|
||||
!et
|
||||
Find these parameters for your data set and compare the fitted function $y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i$ and find the error in the parameters $\beta_0$ and $\beta_1$ as well.
|
||||
How does your parametrization compare with the input data $y_i$? Give a critical discussion of the method.
|
||||
|
||||
Develop now a program where you extract the parameters $\beta_i$ and their corresponding errors using the general Least Square method (LSM), with polynomials of degree $p=2$ and $p=3$. Which polyomials gives the best fit to your data? Here you should prepare an input file with the number of data points in your training data set, the corresponding standard deviation $\sigma_i$ as well as the degree $p$ of the polynomial you want to fit.
|
||||
|
||||
The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients $\beta_i$. A better approach is to use the Singular Value Decomposition (SVD) method discussed in the lecture notes on regression analysis. Here you should use the "SVD programs":"https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/SVD" and use the same input file you used for the general LSM. In the main functions you should specifiy which function you wish to approximate the date with. The example programs discussed above show you examples on how to use these codes as well as how to write the functional form of the polynomial or fitting function.
|
||||
|
||||
|
||||
=== Part d): Introducing Bayesian statistics ===
|
||||
More text to come
|
||||
|
||||
=== Part e): Studying the Ising model or the VMC results with Neural networks ===
|
||||
More text to come
|
||||
|
||||
===== Background literature =====
|
||||
|
||||
If you wish to read more about the Ising model and statistical physics here are three suggestions.
|
||||
|
||||
* "M. Plischke and B. Bergersen":"http://www.worldscientific.com/worldscibooks/10.1142/5660", *Equilibrium Statistical Physics*, World Scientific, see chapters 5 and 6.
|
||||
|
||||
* "D. P. Landau and K. Binder":"http://www.cambridge.org/no/academic/subjects/physics/computational-science-and-modelling/guide-monte-carlo-simulations-statistical-physics-4th-edition?format=HB", *A Guide to Monte Carlo Simulations in Statistical Physics*, Cambridge, see chapters 2,3 and 4.
|
||||
|
||||
* "M. E. J. Newman and T. Barkema":"https://global.oup.com/academic/product/monte-carlo-methods-in-statistical-physics-9780198517979?cc=no&lang=en&", *Monte Carlo Methods in Statistical Physics*, Oxford, see chapters 3 and 4.
|
||||
|
||||
|
||||
===== Introduction to numerical projects =====
|
||||
|
||||
Here follows a brief recipe and recommendation on how to write a report for each
|
||||
project.
|
||||
|
||||
* Give a short description of the nature of the problem and the eventual numerical methods you have used.
|
||||
|
||||
* Describe the algorithm you have used and/or developed. Here you may find it convenient to use pseudocoding. In many cases you can describe the algorithm in the program itself.
|
||||
|
||||
* Include the source code of your program. Comment your program properly.
|
||||
|
||||
* If possible, try to find analytic solutions, or known limits in order to test your program when developing the code.
|
||||
|
||||
* Include your results either in figure form or in a table. Remember to label your results. All tables and figures should have relevant captions and labels on the axes.
|
||||
|
||||
* Try to evaluate the reliabilty and numerical stability/precision of your results. If possible, include a qualitative and/or quantitative discussion of the numerical stability, eventual loss of precision etc.
|
||||
|
||||
* Try to give an interpretation of you results in your answers to the problems.
|
||||
|
||||
* Critique: if possible include your comments and reflections about the exercise, whether you felt you learnt something, ideas for improvements and other thoughts you've made when solving the exercise. We wish to keep this course at the interactive level and your comments can help us improve it.
|
||||
|
||||
* Try to establish a practice where you log your work at the computerlab. You may find such a logbook very handy at later stages in your work, especially when you don't properly remember what a previous test version of your program did. Here you could also record the time spent on solving the exercise, various algorithms you may have tested or other topics which you feel worthy of mentioning.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -267,6 +267,79 @@ meaning that the solution for $\hat{\beta}$ is the one which minimizes the resid
|
||||
|
||||
!eblock
|
||||
|
||||
!split
|
||||
===== Simple regression model =====
|
||||
We are now ready to write our first program which aims at solving the above linear regression equations. We start with data we have produced ourselves, in this case normally distributed random numbers along the $x$-axis. These numbers define then the value of a function $y(x)=4+3x+N(0,1)$. Thereafter we order the $x$ values and employ our linear regression algorithm to set up the best fit.
|
||||
|
||||
!bc pycod
|
||||
# Importing various packages
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
x = 2*np.random.rand(100,1)
|
||||
y = 4+3*x+np.random.randn(100,1)
|
||||
|
||||
xb = np.c_[np.ones((100,1)), x]
|
||||
theta = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
|
||||
xnew = np.array([[0],[2]])
|
||||
xbnew = np.c_[np.ones((2,1)), xnew]
|
||||
ypredict = xbnew.dot(theta)
|
||||
|
||||
plt.plot(xnew, ypredict, "r-")
|
||||
plt.plot(x, y ,'ro')
|
||||
plt.axis([0,2.0,0, 15.0])
|
||||
plt.xlabel(r'$x$')
|
||||
plt.ylabel(r'$y$')
|
||||
plt.title(r'Linear Regression')
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Simple regression model, now using _scikit-learn_ =====
|
||||
We can repeat the above algorithm using _scikit-learn_ as follows.
|
||||
!bc pycod
|
||||
# Importing various packages
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.linear_model import LinearRegression
|
||||
|
||||
x = 2*np.random.rand(100,1)
|
||||
y = 4+3*x+np.random.randn(100,1)
|
||||
linreg = LinearRegression()
|
||||
linreg.fit(x,y)
|
||||
xnew = np.array([[0],[2]])
|
||||
ypredict = linreg.predict(xnew)
|
||||
|
||||
plt.plot(xnew, ypredict, "r-")
|
||||
plt.plot(x, y ,'ro')
|
||||
plt.axis([0,2.0,0, 15.0])
|
||||
plt.xlabel(r'$x$')
|
||||
plt.ylabel(r'$y$')
|
||||
plt.title(r'Random numbers ')
|
||||
plt.show()
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Correlations and the quality of our results =====
|
||||
In order to test the quality of our fit, there are several measures which can be implemented. One is the so-called
|
||||
correlation function defined as
|
||||
!bt
|
||||
\[
|
||||
\mathrm{Corr}(X,Y) = \frac{\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum_{i=1}^n(x_i-\overline{x})^2}\sqrt{\sum_{i=1}^n(y_i-\overline{y})^2} }
|
||||
\]
|
||||
!et
|
||||
|
||||
Another quantity is the autocorrelation function we discussed in our chapter on statistical analysis.
|
||||
|
||||
!split
|
||||
===== The $\chi^2$ function =====
|
||||
!bblock
|
||||
@@ -433,6 +506,83 @@ The LSM suffers often from both being underdetermined and overdetermined in the
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Simple regression model with gradient descent =====
|
||||
Add info about the equations, play around with different learning rates
|
||||
!bc pycod
|
||||
# Importing various packages
|
||||
from math import exp, sqrt
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
x = 2*np.random.rand(100,1)
|
||||
y = 4+3*x+np.random.randn(100,1)
|
||||
|
||||
xb = np.c_[np.ones((100,1)), x]
|
||||
theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
|
||||
print(theta_linreg)
|
||||
theta = np.random.randn(2,1)
|
||||
|
||||
eta = 0.1
|
||||
Niterations = 1000
|
||||
m = 100
|
||||
|
||||
for iter in range(Niterations):
|
||||
gradients = 2.0/m*xb.T.dot(xb.dot(theta)-y)
|
||||
theta -= eta*gradients
|
||||
|
||||
print(theta)
|
||||
xnew = np.array([[0],[2]])
|
||||
xbnew = np.c_[np.ones((2,1)), xnew]
|
||||
ypredict = xbnew.dot(theta)
|
||||
ypredict2 = xbnew.dot(theta_linreg)
|
||||
plt.plot(xnew, ypredict, "r-")
|
||||
plt.plot(xnew, ypredict2, "b-")
|
||||
plt.plot(x, y ,'ro')
|
||||
plt.axis([0,2.0,0, 15.0])
|
||||
plt.xlabel(r'$x$')
|
||||
plt.ylabel(r'$y$')
|
||||
plt.title(r'Random numbers ')
|
||||
plt.show()
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Simple regression model with stochastic gradient descent =====
|
||||
Add info about the equations, play around with different learning rates
|
||||
!bc pycod
|
||||
# Importing various packages
|
||||
from math import exp, sqrt
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.linear_model import SGDRegressor
|
||||
|
||||
x = 2*np.random.rand(100,1)
|
||||
y = 4+3*x+np.random.randn(100,1)
|
||||
|
||||
xb = np.c_[np.ones((100,1)), x]
|
||||
theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
|
||||
print(theta_linreg)
|
||||
sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
|
||||
sgdreg.fit(x,y.ravel())
|
||||
print(sgdreg.intercept_, sgdreg.coef_)
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== The singular value decompostion =====
|
||||
!bblock
|
||||
|
||||
+25
-25
@@ -44,34 +44,34 @@ system doconce split_html $html.html --method=space10
|
||||
# Bootstrap style
|
||||
html=${name}-bs
|
||||
system doconce format html $name --html_style=bootstrap --pygments_html_style=default --html_admon=bootstrap_panel --html_output=$html $opt
|
||||
system doconce split_html $html.html --method=split --pagination --nav_button=bottom
|
||||
#system doconce split_html $html.html --method=split --pagination --nav_button=bottom
|
||||
|
||||
# IPython notebook
|
||||
#system doconce format ipynb $name $opt
|
||||
system doconce format ipynb $name $opt
|
||||
|
||||
#--# LaTeX Beamer slides
|
||||
#--beamertheme=red_plain
|
||||
#--system doconce format pdflatex $name --latex_title_layout=beamer --latex_table_format=footnotesize $opt
|
||||
#--system doconce ptex2tex $name envir=minted
|
||||
#--# Add special packages
|
||||
#--doconce subst "% Add user's preamble" "\g<1>\n\\usepackage{simplewick}" $name.tex
|
||||
#--system doconce slides_beamer $name --beamer_slide_theme=$beamertheme
|
||||
#--system pdflatex -shell-escape ${name}
|
||||
#--system pdflatex -shell-escape ${name}
|
||||
#--cp $name.pdf ${name}-beamer.pdf
|
||||
#--cp $name.tex ${name}-beamer.tex
|
||||
#--
|
||||
#--# Handouts
|
||||
#--system doconce format pdflatex $name --latex_title_layout=beamer --latex_table_format=footnotesize $opt
|
||||
#--system doconce ptex2tex $name envir=minted
|
||||
#--# Add special packages
|
||||
#--doconce subst "% Add user's preamble" "\g<1>\n\\usepackage{simplewick}" $name.tex
|
||||
#--system doconce slides_beamer $name --beamer_slide_theme=red_shadow --handout
|
||||
#--system pdflatex -shell-escape $name
|
||||
#--pdflatex -shell-escape $name
|
||||
#--pdflatex -shell-escape $name
|
||||
#--pdfnup --nup 2x3 --frame true --delta "1cm 1cm" --scale 0.9 --outfile ${name}-beamer-handouts2x3.pdf ${name}.pdf
|
||||
#--rm -f ${name}.pdf
|
||||
# LaTeX Beamer slides
|
||||
beamertheme=red_plain
|
||||
system doconce format pdflatex $name --latex_title_layout=beamer --latex_table_format=footnotesize $opt
|
||||
system doconce ptex2tex $name envir=minted
|
||||
# Add special packages
|
||||
doconce subst "% Add user's preamble" "\g<1>\n\\usepackage{simplewick}" $name.tex
|
||||
system doconce slides_beamer $name --beamer_slide_theme=$beamertheme
|
||||
system pdflatex -shell-escape ${name}
|
||||
system pdflatex -shell-escape ${name}
|
||||
cp $name.pdf ${name}-beamer.pdf
|
||||
cp $name.tex ${name}-beamer.tex
|
||||
|
||||
# Handouts
|
||||
system doconce format pdflatex $name --latex_title_layout=beamer --latex_table_format=footnotesize $opt
|
||||
system doconce ptex2tex $name envir=minted
|
||||
# Add special packages
|
||||
doconce subst "% Add user's preamble" "\g<1>\n\\usepackage{simplewick}" $name.tex
|
||||
system doconce slides_beamer $name --beamer_slide_theme=red_shadow --handout
|
||||
system pdflatex -shell-escape $name
|
||||
pdflatex -shell-escape $name
|
||||
pdflatex -shell-escape $name
|
||||
pdfnup --nup 2x3 --frame true --delta "1cm 1cm" --scale 0.9 --outfile ${name}-beamer-handouts2x3.pdf ${name}.pdf
|
||||
rm -f ${name}.pdf
|
||||
|
||||
# Ordinary plain LaTeX document
|
||||
rm -f *.aux # important after beamer
|
||||
|
||||
Reference in New Issue
Block a user