FYS-STK4155/doc/pub/BM/ipynb/BM.ipynb

Machine Learning and Boltzmann machines with applications

Morten Hjorth-Jensen, Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University and Department of Physics, University of Oslo, Norway

Date: Nov 29, 2018

Copyright 1999-2018, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license

Types of Machine Learning, a repetition

The approaches to machine learning are many, but are often split into two main categories. 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, solely from rewards and punishment.

Another way to categorize machine learning tasks is to consider the desired output of a system. Some of the most common tasks are:

• Classification: Outputs are divided into two or more classes. The goal is to produce a model that assigns inputs into one of these classes. An example is to identify digits based on pictures of hand-written ones. Classification is typically supervised learning.

• Regression: Finding a functional relationship between an input data set and a reference data set. The goal is to construct a function that maps input data to continuous output values.

• Clustering: Data are divided into groups with certain common traits, without knowing the different groups beforehand. It is thus a form of unsupervised learning.

• Other unsupervised learning algortihms, here Boltzmann machines

Why Boltzmann machines?

What is known as restricted Boltzmann Machines (RMB) have received a lot of attention lately. One of the major reasons is that they can be stacked layer-wise to build deep neural networks that capture complicated statistics.

The original RBMs had just one visible layer and a hidden layer, but recently so-called Gaussian-binary RBMs have gained quite some popularity in imaging since they are capable of modeling continuous data that are common to natural images.

Furthermore, they have been used to solve complicated quantum mechanical many-particle problems or classical statistical physics problems like the Ising and Potts classes of models.

An intermediate step, the Hopfield network and links to the Ising and Potts models

More material on Hopfield networks will come here later.

A brief review on Markov Chains, Metropolis and Gibbs sampling

• We want to study a physical system which evolves towards equilibrium, from given initial conditions.

• We start with a PDF w(x_0,t_0) and we want to understand how the system evolves with time.

• We want to reach a situation where after a given number of time steps we obtain a steady state. This means that the system reaches its most likely state (equilibrium situation)

• Our PDF is normally a multidimensional object whose normalization constant is impossible to find.

• Analytical calculations from w(x,t) are not possible.

• To sample directly from from w(x,t) is not possible/difficult.

• The transition probability W is also not known.

• How can we establish that we have reached a steady state? Sounds impossible!

Use Markov chain Monte Carlo

Brownian motion and Markov processes

A Markov process is a random walk with a selected probability for making a move. The new move is independent of the previous history of the system.

The Markov process is used repeatedly in Monte Carlo simulations in order to generate new random states.

The reason for choosing a Markov process is that when it is run for a long enough time starting with a random state, we will eventually reach the most likely state of the system.

In thermodynamics, this means that after a certain number of Markov processes we reach an equilibrium distribution.

This mimicks the way a real system reaches its most likely state at a given temperature of the surroundings.

Brownian motion and Markov processes, Ergodicity and Detailed balance

To reach this distribution, the Markov process needs to obey two important conditions, that of ergodicity and detailed balance. These conditions impose then constraints on our algorithms for accepting or rejecting new random states.

The Metropolis algorithm discussed here abides to both these constraints.

The Metropolis algorithm is widely used in Monte Carlo simulations and the understanding of it rests within the interpretation of random walks and Markov processes.

Brownian motion and Markov processes, jargon

In a random walk one defines a mathematical entity called a walker, whose attributes completely define the state of the system in question.

The state of the system can refer to any physical quantities, from the vibrational state of a molecule specified by a set of quantum numbers, to the brands of coffee in your favourite supermarket.

The walker moves in an appropriate state space by a combination of deterministic and random displacements from its previous position.

This sequence of steps forms a chain.

Brownian motion and Markov processes, sequence of ingredients

• We want to study a physical system which evolves towards equilibrium, from given initial conditions.

• Markov chains are intimately linked with the physical process of diffusion.

• From a Markov chain we can then derive the conditions for detailed balance and ergodicity. These are the conditions needed for obtaining a steady state.

• The widely used algorithm for doing this is the so-called Metropolis algorithm, in its refined form the Metropolis-Hastings algorithm.

Applications: almost every field in science

• Financial engineering, see for example Patriarca et al, Physica 340, page 334 (2004).

• Neuroscience, see for example Lipinski, Physics Medical Biology 35, page 441 (1990) or Farnell and Gibson, Journal of Computational Physics 208, page 253 (2005)

• Tons of applications in physics

• and chemistry

• and biology, medicine

• Nobel prize in economy to Black and Scholes

\frac{\partial V}{\partial t}+\frac{1}{2}\sigma^{2}S^{2}\frac{\partial^{2} V}{\partial S^{2}}+rS\frac{\partial V}{\partial S}-rV=0.

The Black and Scholes equation is a partial differential equation, which describes the price of the option over time. It is a diffusion equation with a random term.

The list of applications is endless.

Markov processes

A Markov process allows in principle for a microscopic description of Brownian motion. As with the random walk studied in the previous section, we consider a particle which moves along the $x$-axis in the form of a series of jumps with step length \Delta x = l. Time and space are discretized and the subsequent moves are statistically independent, i.e., the new move depends only on the previous step and not on the results from earlier trials. We start at a position x=jl=j\Delta x and move to a new position x =i\Delta x during a step \Delta t=\epsilon, where i\ge 0 and j\ge 0 are integers. The original probability distribution function (PDF) of the particles is given by
w_i(t=0) where i refers to a specific position on the grid in

The function w_i(t=0) is now the discretized version of w(x,t). We can regard the discretized PDF as a vector.

Markov processes

For the Markov process we have a transition probability from a position x=jl to a position x=il given by

W_{ij}(\epsilon)=W(il-jl,\epsilon)=\left\{\begin{array}{cc}\frac{1}{2} & |i-j| = 1\\
             0 & \mathrm{else} \end{array} \right. ,

where W_{ij} is normally called the transition probability and we can represent it, see below, as a matrix. Here we have specialized to a case where the transition probability is known.

Our new PDF w_i(t=\epsilon) is now related to the PDF at t=0 through the relation

w_i(t=\epsilon) =\sum_{j} W(j\rightarrow i)w_j(t=0).

This equation represents the discretized time-development of an original PDF with equal probability of jumping left or right.

Markov processes, the probabilities

Since both W and w represent probabilities, they have to be normalized, i.e., we require that at each time step we have

\sum_i w_i(t) = 1,

and

\sum_j W(j\rightarrow i) = 1,

which applies for all $j$-values. The further constraints are 0 \le W_{ij} \le 1 and 0 \le w_{j} \le 1. Note that the probability for remaining at the same place is in general not necessarily equal zero.

Markov processes

The time development of our initial PDF can now be represented through the action of the transition probability matrix applied n times. At a time t_n=n\epsilon our initial distribution has developed into

w_i(t_n) = \sum_jW_{ij}(t_n)w_j(0),

and defining

W(il-jl,n\epsilon)=(W^n(\epsilon))_{ij}

we obtain

w_i(n\epsilon) = \sum_j(W^n(\epsilon))_{ij}w_j(0),

or in matrix form

\begin{equation} \label{eq:wfinal} \tag{1}
   \hat{w}(n\epsilon) = \hat{W}^n(\epsilon)\hat{w}(0).
\end{equation}

An Illustrative Example

The following simple example may help in understanding the meaning of the transition matrix \hat{W} and the vector \hat{w}. Consider the 4\times 4 matrix \hat{W}

\hat{W} = \left(\begin{array}{cccc} 1/4 & 1/9 & 3/8 & 1/3 \\                   
                                       2/4 & 2/9 & 0 & 1/3\\                   
                                       0   & 1/9 & 3/8 & 0\\
                                       1/4 & 5/9&  2/8 & 1/3 \end{array} \right),

and we choose our initial state as

\hat{w}(t=0)=  \left(\begin{array}{c} 1\\                   
                                 0\\
                                 0 \\                   
                                 0 \end{array} \right).

An Illustrative Example

We note that both the vector and the matrix are properly normalized. Summing the vector elements gives one and summing over columns for the matrix results also in one. Furthermore, the largest eigenvalue is one. We act then on \hat{w} with \hat{W}. The first iteration is

\hat{w}(t=\epsilon) = \hat{W}\hat{w}(t=0),

resulting in

\hat{w}(t=\epsilon)=  \left(\begin{array}{c} 1/4\\                   
                                1/2 \\
                                0 \\                   
                                1/4 \end{array} \right).

An Illustrative Example, next step

The next iteration results in

\hat{w}(t=2\epsilon) = \hat{W}\hat{w}(t=\epsilon),

resulting in

\hat{w}(t=2\epsilon)=  \left(\begin{array}{c} 0.201389\\
   0.319444 \\
   0.055556 \\
   0.423611 \end{array} \right).

Note that the vector \hat{w} is always normalized to 1.

An Illustrative Example, the steady state

We find the steady state of the system by solving the set of equations

w(t=\infty) = Ww(t=\infty),

which is an eigenvalue problem with eigenvalue equal to one! This set of equations reads

W_{11}w_1(t=\infty) +W_{12}w_2(t=\infty) +W_{13}w_3(t=\infty)+ W_{14}w_4(t=\infty)=w_1(t=\infty) \nonumber

W_{21}w_1(t=\infty) + W_{22}w_2(t=\infty) + W_{23}w_3(t=\infty)+ W_{24}w_4(t=\infty)=w_2(t=\infty) \nonumber

W_{31}w_1(t=\infty) + W_{32}w_2(t=\infty) + W_{33}w_3(t=\infty)+ W_{34}w_4(t=\infty)=w_3(t=\infty) \nonumber

W_{41}w_1(t=\infty) + W_{42}w_2(t=\infty) + W_{43}w_3(t=\infty)+ W_{44}w_4(t=\infty)=w_4(t=\infty) \nonumber

\begin{equation} 
\label{_auto1} \tag{2}
\end{equation}

with the constraint that

\sum_i w_i(t=\infty) = 1,

yielding as solution

\hat{w}(t=\infty)=  \left(\begin{array}{c}0.244318 \\                   
                                 0.319602 \\  0.056818 \\  0.379261 \end{array} \right).

An Illustrative Example, iterative steps

The table here demonstrates the convergence as a function of the number of iterations or time steps. After twelve iterations we have reached the exact value with six leading digits.

Iteration	$w_1$	$w_2$	$w_3$	$w_4$
0	1.000000	0.000000	0.000000	0.000000
1	0.250000	0.500000	0.000000	0.250000
2	0.201389	0.319444	0.055556	0.423611
3	0.247878	0.312886	0.056327	0.382909
4	0.245494	0.321106	0.055888	0.377513
5	0.243847	0.319941	0.056636	0.379575
6	0.244274	0.319547	0.056788	0.379391
7	0.244333	0.319611	0.056801	0.379255
8	0.244314	0.319610	0.056813	0.379264
9	0.244317	0.319603	0.056817	0.379264
10	0.244318	0.319602	0.056818	0.379262
11	0.244318	0.319602	0.056818	0.379261
12	0.244318	0.319602	0.056818	0.379261
$\hat{w}(t=\infty)$	0.244318	0.319602	0.056818	0.379261

An Illustrative Example, what does it mean?

We have after $t$-steps

\hat{w}(t) = \hat{W}^t\hat{w}(0),

with \hat{w}(0) the distribution at t=0 and \hat{W} representing the transition probability matrix.

An Illustrative Example, understanding the basics

We can always expand \hat{w}(0) in terms of the right eigenvectors \hat{v} of \hat{W} as

\hat{w}(0)  = \sum_i\alpha_i\hat{v}_i,

resulting in

\hat{w}(t) = \hat{W}^t\hat{w}(0)=\hat{W}^t\sum_i\alpha_i\hat{v}_i=
\sum_i\lambda_i^t\alpha_i\hat{v}_i,

with \lambda_i the i^{\mathrm{th}} eigenvalue corresponding to
the eigenvector \hat{v}_i.

If we assume that \lambda_0 is the largest eigenvector we see that in the limit t\rightarrow \infty, \hat{w}(t) becomes proportional to the corresponding eigenvector \hat{v}_0. This is our steady state or final distribution.

The Metropolis Algorithm and Detailed Balance

Let us recapitulate some of our results about Markov chains and random walks.

• The time development of our PDF w(t), after

one time-step from t=0 is given by

w_i(t=\epsilon) = W(j\rightarrow i)w_j(t=0).

This equation represents the discretized time-development of an original PDF. We can rewrite this as a

w_i(t=\epsilon) = W_{ij}w_j(t=0).

with the transition matrix W for a random walk given by

W_{ij}(\epsilon)=W(il-jl,\epsilon)=\left\{\begin{array}{cc}\frac{1}{2} & |i-j| = 1\\
                                             0 & \mathrm{else} \end{array} \right.

The Metropolis Algorithm and Detailed Balance

We call W_{ij} for the transition probability and we represent it as a matrix.

• Both W and w represent probabilities and they have to be normalized, meaning that at each time step we have

\sum_i w_i(t) = 1,

and

\sum_j W(j\rightarrow i) = 1.

Here we have written the previous matrix W_{ij}=W(j\rightarrow i).

The Metropolis Algorithm and Detailed Balance

The further constraints are 0 \le W_{ij} \le 1 and 0 \le w_{j} \le 1.

• We can thus write the action of W as

w_i(t+1) = \sum_jW_{ij}w_j(t),

or as vector-matrix relation

\hat{w}(t+1) = \hat{W\hat{w}}(t),

and if we have that ||\hat{w}(t+1)-\hat{w}(t)||\rightarrow 0, we say that we have reached the most likely state of the system, the so-called steady state or equilibrium state.

The Metropolis Algorithm and Detailed Balance

Another way of phrasing this is

\begin{equation} 
w(t=\infty) = Ww(t=\infty). 
\label{_auto2} \tag{3}
\end{equation}

The Metropolis Algorithm and Detailed Balance

The question then is how can we model anything under such a severe lack of knowledge? The Metropolis algorithm comes to our rescue here. Since W(j\rightarrow i) is unknown, we model it as the product of two probabilities, a probability for accepting the proposed move from the state j to the state j, and a probability for making the transition to the state i being in the state j. We label these probabilities A(j\rightarrow i) and T(j\rightarrow i), respectively. Our total transition probability is then

W(j\rightarrow i)=T(j\rightarrow i)A(j\rightarrow i).

The algorithm can then be expressed as

• We make a suggested move to the new state i with some transition or moving probability T_{j\rightarrow i}.

• We accept this move to the new state with an acceptance probability A_{j \rightarrow i}. The new state i is in turn used as our new starting point for the next move. We reject this proposed moved with a 1-A_{j\rightarrow i} and the original state j is used again as a sample.

The Metropolis Algorithm and Detailed Balance

We wish to derive the required properties of the probabilities T and A such that w_i^{(t\rightarrow \infty)} \rightarrow w_i, starting from any distribution, will lead us to the correct distribution.

We can now derive the dynamical process towards equilibrium. To obtain this equation we note that after t time steps the probability for being in a state i is related to the probability of being in a state j and performing a transition to the new state together with the probability of actually being in the state i and making a move to any of the possible states j from the previous time step.

The Metropolis Algorithm and Detailed Balance

We can express this as, assuming that T and A are time-independent,

w_i(t+1) = \sum_j \left [
w_j(t)T_{j\rightarrow i} A_{j\rightarrow i} 
+w_i(t)T_{i\rightarrow j}\left ( 1- A_{i\rightarrow j} \right)
\right ] \,.

The Metropolis Algorithm and Detailed Balance

All probabilities are normalized, meaning that \sum_j T_{i\rightarrow j} = 1. Using the latter, we can rewrite the previous equation as

w_i(t+1) = w_i(t) +
 \sum_j \left [
w_j(t)T_{j\rightarrow i} A_{j\rightarrow i} 
-w_i(t)T_{i\rightarrow j}A_{i\rightarrow j}\right ] \,,

which can be rewritten as

w_i(t+1)-w_i(t) =  \sum_j \left [w_j(t)T_{j\rightarrow i} A_{j\rightarrow i} 
-w_i(t)T_{i\rightarrow j}A_{i\rightarrow j}\right ] .

The Metropolis Algorithm and Detailed Balance

The last equation is very similar to the so-called Master equation, which relates the temporal dependence of a PDF w_i(t) to various transition rates. The equation can be derived from the so-called Chapman-Einstein-Enskog-Kolmogorov equation. The equation is given as

\begin{equation}
\label{eq:masterequation} \tag{4}
\frac{d w_i(t)}{dt} = \sum_j\left[ W(j\rightarrow i)w_j-W(i\rightarrow j)w_i\right],
\end{equation}

which simply states that the rate at which the systems moves from a state j to a final state i (the first term on the right-hand side of the last equation) is balanced by the rate at which the system undergoes transitions from the state i to a state j (the second term). If we have reached the so-called steady state, then the temporal development is zero. This means that in equilibrium we have

\frac{d w_i(t)}{dt} = 0.

The Metropolis Algorithm and Detailed Balance

In the limit t\rightarrow \infty we require that the two distributions w_i(t+1)=w_i and w_i(t)=w_i and we have

\sum_j w_jT_{j\rightarrow i} A_{j\rightarrow i}= \sum_j w_iT_{i\rightarrow j}A_{i\rightarrow j},

which is the condition for balance when the most likely state (or steady state) has been reached. We see also that the right-hand side can be rewritten as

\sum_j w_iT_{i\rightarrow j}A_{i\rightarrow j}= \sum_j w_iW_{i\rightarrow j},

and using the property that \sum_j W_{i\rightarrow j}=1, we can rewrite our equation as

w_i= \sum_j w_jT_{j\rightarrow i} A_{j\rightarrow i}= \sum_j w_j W_{j\rightarrow i},

which is nothing but the standard equation for a Markov chain when the steady state has been reached.

The Metropolis Algorithm and Detailed Balance

However, the condition that the rates should equal each other is in general not sufficient to guarantee that we, after many simulations, generate the correct distribution. We may risk to end up with so-called cyclic solutions. To avoid this we therefore introduce an additional condition, namely that of detailed balance

W(j\rightarrow i)w_j= W(i\rightarrow j)w_i.

These equations were derived by Lars Onsager when studying irreversible processes. At equilibrium detailed balance gives thus

\frac{W(j\rightarrow i)}{W(i\rightarrow j)}=\frac{w_i}{w_j}.

Rewriting the last equation in terms of our transition probabilities T and acceptance probobalities A we obtain

w_j(t)T_{j\rightarrow i}A_{j\rightarrow i}= w_i(t)T_{i\rightarrow j}A_{i\rightarrow j}.

The Metropolis Algorithm and Detailed Balance

Since we normally have an expression for the probability distribution functions w_i, we can rewrite the last equation as

\frac{T_{j\rightarrow i}A_{j\rightarrow i}}{T_{i\rightarrow j}A_{i\rightarrow j}}= \frac{w_i}{w_j}.

The Metropolis Algorithm and Detailed Balance

In statistical physics this condition ensures that it is e.g., the Boltzmann distribution which is generated when equilibrium is reached.

We introduce now the Boltzmann distribution

w_i= \frac{\exp{(-\beta(E_i))}}{Z},

which states that the probability of finding the system in a state i with energy E_i at an inverse temperature \beta = 1/k_BT is w_i\propto \exp{(-\beta(E_i))}. The denominator Z is a normalization constant which ensures that the sum of all probabilities is normalized to one. It is defined as the sum of probabilities over all microstates j of the system

Z=\sum_j \exp{(-\beta(E_i))}.

The Metropolis Algorithm and Detailed Balance

From the partition function we can in principle generate all interesting quantities for a given system in equilibrium with its surroundings at a temperature T.

With the probability distribution given by the Boltzmann distribution we are now in a position where we can generate expectation values for a given variable A through the definition

\langle A \rangle = \sum_jA_jw_j=
    \frac{\sum_jA_j\exp{(-\beta(E_j)}}{Z}.

In general, most systems have an infinity of microstates making thereby the computation of Z practically impossible and a brute force Monte Carlo calculation over a given number of randomly selected microstates may therefore not yield those microstates which are important at equilibrium. To select the most important contributions we need to
use the condition for detailed balance. Since this is just given by the ratios of probabilities, we never need to evaluate the partition function Z.

The Metropolis Algorithm and Detailed Balance

For the Boltzmann distribution, detailed balance results in

\frac{w_i}{w_j}= \exp{(-\beta(E_i-E_j))}.

Let us now specialize to a system whose energy is defined by the orientation of single spins. Consider the state i, with given energy E_i represented by the following N spins

\begin{array}{cccccccccc}
\uparrow&\uparrow&\uparrow&\dots&\uparrow&\downarrow&\uparrow&\dots&\uparrow&\downarrow\\
1&2&3&\dots& k-1&k&k+1&\dots&N-1&N\end{array}

Warning:

Output truncated. This notebook contains too many cells to display efficiently.