From 3d7ab1cd12809390e70dcccc70c7a9a63ea614f5 Mon Sep 17 00:00:00 2001 From: mhjensen Date: Mon, 19 Nov 2018 22:23:13 +0100 Subject: [PATCH] decision tree not added properly --- doc/pub/BM/html/BM-bs.html | 1069 +------------ doc/pub/BM/html/BM-reveal.html | 947 +----------- doc/pub/BM/html/BM-solarized.html | 960 +----------- doc/pub/BM/html/BM.html | 960 +----------- doc/pub/BM/ipynb/BM.ipynb | 1005 +----------- doc/pub/BM/ipynb/ipynb-BM-src.tar.gz | Bin 429509 -> 70532 bytes doc/pub/BM/pdf/BM-beamer-handouts2x3.pdf | Bin 1176339 -> 612544 bytes doc/pub/BM/pdf/BM-beamer.pdf | Bin 1158940 -> 594066 bytes doc/pub/BM/pdf/BM-minted.pdf | Bin 1161785 -> 632243 bytes doc/pub/DecisionTrees/html/data/zoo.csv | 101 ++ .../reveal.js/css/theme/template/mixins.scss | 29 + .../css/theme/template/settings.scss | 34 + .../reveal.js/css/theme/template/theme.scss | 171 +++ .../DecisionTrees/html/reveal.js/demo.html | 410 +++++ .../reveal.js/plugin/multiplex/package.json | 19 + .../html/reveal.js/test/simple.md | 12 + .../test/test-markdown-external.html | 36 + .../reveal.js/test/test-markdown-external.js | 24 + .../reveal.js/test/test-markdown-options.html | 41 + .../reveal.js/test/test-markdown-options.js | 26 + .../DecisionTrees-checkpoint.ipynb | 416 +++++ .../DecisionTrees/ipynb/DecisionTrees.ipynb | 1365 +++++++++++++++++ doc/pub/DecisionTrees/ipynb/data/zoo.csv | 101 ++ doc/src/BoltzmannMachines/BM.do.txt | 652 +------- 24 files changed, 3080 insertions(+), 5298 deletions(-) create mode 100644 doc/pub/DecisionTrees/html/data/zoo.csv create mode 100644 doc/pub/DecisionTrees/html/reveal.js/css/theme/template/mixins.scss create mode 100644 doc/pub/DecisionTrees/html/reveal.js/css/theme/template/settings.scss create mode 100644 doc/pub/DecisionTrees/html/reveal.js/css/theme/template/theme.scss create mode 100644 doc/pub/DecisionTrees/html/reveal.js/demo.html create mode 100644 doc/pub/DecisionTrees/html/reveal.js/plugin/multiplex/package.json create mode 100644 doc/pub/DecisionTrees/html/reveal.js/test/simple.md create mode 100644 doc/pub/DecisionTrees/html/reveal.js/test/test-markdown-external.html create mode 100644 doc/pub/DecisionTrees/html/reveal.js/test/test-markdown-external.js create mode 100644 doc/pub/DecisionTrees/html/reveal.js/test/test-markdown-options.html create mode 100644 doc/pub/DecisionTrees/html/reveal.js/test/test-markdown-options.js create mode 100644 doc/pub/DecisionTrees/ipynb/.ipynb_checkpoints/DecisionTrees-checkpoint.ipynb create mode 100644 doc/pub/DecisionTrees/ipynb/DecisionTrees.ipynb create mode 100644 doc/pub/DecisionTrees/ipynb/data/zoo.csv diff --git a/doc/pub/BM/html/BM-bs.html b/doc/pub/BM/html/BM-bs.html index a2a6f79f8..a4996e476 100644 --- a/doc/pub/BM/html/BM-bs.html +++ b/doc/pub/BM/html/BM-bs.html @@ -6,9 +6,9 @@ Automatically generated HTML file from DocOnce source - + -Solving quantum mechanical problems with Machine Learning and Boltzmann machines +Machine Learning and Boltzmann machines @@ -40,95 +40,44 @@ Automatically generated HTML file from DocOnce source @@ -158,7 +107,7 @@ MathJax.Hub.Config({ - Solving quantum mechanical problems with Machine Learning and Boltzmann machines + Machine Learning and Boltzmann machines
@@ -166,63 +115,32 @@ MathJax.Hub.Config({
  • Contents
    • @@ -240,7 +158,7 @@ MathJax.Hub.Config({
    -

    Solving quantum mechanical problems with Machine Learning and Boltzmann machines

    +

    Machine Learning and Boltzmann machines

    @@ -255,7 +173,7 @@ MathJax.Hub.Config({

    Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University and Department of Physics, University of Oslo, Norway

    -

    Jul 26, 2018

    +

    Nov 19, 2018


    @@ -263,18 +181,10 @@ MathJax.Hub.Config({ -

    What is this talk about?

    +

    Unsupervised learning, ovrarching aims

    -The main aim is to give you a short and pedestrian introduction to how we can use Machine Learning methods -to solve quantum mechanical many-body problems. And why this could be of interest. - -

    -The hope is that after this talk you have gotten the basic ideas to get you started. Peeping into https://github.com/mhjensenseminars/MachineLearningTalk, you'll find a Jupyter notebook, slides, codes etc that will allow you to reproduce the simulations discussed here, and perhaps run your own very first calculations. - -

    -Furthermore, since we will also talk about Monte Carlo simulations, the exciting idea of Monte Carlo optimization and Machine Learning can be used to mesh strategies for reducing the curse of dimensionality.

    @@ -284,196 +194,7 @@ Furthermore, since we will also talk about Monte Carlo simulations, the exciting

    -

    Why?

    - -

    -How can we avoid the dimensionality curse? Many possibilities - -

      -
    1. smarter basis functions
      2. -
    2. resummation of specific correlations
      3. -
    3. stochastic sampling of high-lying states (stochastic FCI, CC and soon SRG)
      4. -
    4. many more
      5. -
    - -Machine Learning and Quantum Computing hold also great promise in tackling the -ever increasing dimensionalities. Here we will focus on Machine Learning. - -

    - - -

    Overview

    -
    -
    -

    - -

      -
    • Short intro to Machine Learning
      • -
    • Variational Monte Carlo (Markov Chain Monte Carlo, \( \mathrm{MC}^2 \)) and two-electron quantum dots, solving quantum mechanical problems in a stochastic way. It will serve as our motivation for switching to Machine Learning.
      • -
    • From Variational Monte Carlo to Boltzmann Machines and Machine Learning
      • -
    -
    -
    - - -

    - - -

    What are the Machine Learning calculations here based on?

    - -

    -This work is inspired by the idea of representing the wave function with -a restricted Boltzmann machine (RBM), presented recently by G. Carleo and M. Troyer, Science 355, Issue 6325, pp. 602-606 (2017). They -named such a wave function/network a neural network quantum state (NQS). In their article they apply it to the quantum mechanical -spin lattice systems of the Ising model and Heisenberg model, with -encouraging results. - -

    -Thanks to Vilde Flugsrud and Alfred Alocias Mariadason for many discussions and interpretations of the above paper. - -

    - - -

    A new world

    -
    -
    -

    -Machine learning (ML) 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 ML are - -

    - -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. -
    -
    - - -

    - - -

    Lots of room for creativity

    -
    -
    -

    -Not all the -algorithms and methods can be given a rigorous mathematical -justification, opening up thereby for experimenting -and trial and error and thereby exciting new developments. -

    -
    - - -

    -

    -
    -

    -A solid command of linear algebra, multivariate theory, -probability theory, statistical data analysis, -understanding errors and Monte Carlo methods is important in order to understand many of the -various algorithms and methods. -

    -
    - - -

    -

    -
    -

    -A familiarity with ML is almost becoming a prerequisite for many of the most exciting employment -opportunities. Add quantum computing and there you are. -

    -
    - - -

    - - -

    Some members of the ML family

    - -

    -

    -
    -

    - -

      -
    1. Linear regression and its variants, in essence polynomial regression
      2. -
    2. Decision tree algorithms, from simpler to more complex ones like random forests
      3. -
    3. Nearest neighbors models
      4. -
    4. Bayesian statistics and regression
      5. -
    5. Support vector machines and finally various variants of
      6. -
    6. Artifical neural networks and deep learning
      7. -
    7. ....
      8. -
    -
    -
    - - -

    - - -

    What are the basic ingredients?

    -
    -
    -

    -Almost every problem in ML and data science starts with the same ingredients: - -

      -
    • The dataset \( \mathbf{x} \) (could be some observable quantity of the system we are studying)
      • -
    • A model which is a function of a set of parameters \( \mathbf{\alpha} \) that relates to the dataset, say a likelihood function \( p(\mathbf{x}\vert \mathbf{\alpha}) \) or just a simple model \( f(\mathbf{\alpha}) \)
      • -
    • A so-called cost function \( \mathcal{C} (\mathbf{x}, f(\mathbf{\alpha})) \) which allows us to decide how well our model represents the dataset.
      • -
    - -We seek to minimize the function \( \mathcal{C} (\mathbf{x}, f(\mathbf{\alpha})) \) by finding the parameter values which minimize \( \mathcal{C} \). This leads to various minimization algorithms. -
    -
    - - -

    - - -

    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 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. - -

    -You will notice however that many of the basic ideas discussed do come from Physics! -

    -
    - - -

    - - -

    Types of Machine Learning

    +

    Types of Machine Learning

    @@ -503,631 +224,7 @@ Some of the most common tasks are:

    -

    References

    -
    -
    -

    - -

    -
    -
    - - -

    - - -

    Another interesting article

    - -

    -Here we will use so-called reduced Boltzmann Machines to simulate quantum many-body problems. For Monte Carlo aficionados, there is a very close similarity with what are called shadow wave functions, see the work of Pederiva and Kalos and collaborators, Phys Rev. E 90, 053304 (2014). - -

    - - -

    Just for fun illustration of scikit-learn: Hudson bay data

    -
    -
    -

    -Most mammalian predators rely on a variety of prey, which complicates mathematical modeling; however, a few predators have become highly specialized and seek almost exclusively a single prey species. An example of this simplified predator-prey interaction is seen in Canadian northern forests, where the populations of the lynx and the snowshoe hare are intertwined in a life and death struggle. - -

    -One reason that this particular system has been so extensively studied is that the Hudson Bay company kept careful records of all furs from the early 1800s into the 1900s. The records for the furs collected by the Hudson Bay company showed distinct oscillations (approximately 12 year periods), suggesting that these species caused almost periodic fluctuations of each other's populations. The table here shows data from 1900 to 1920. - -

    -

    -
    - - -

    - - -

    Hares and lynx in Hudson bay from 1900 to 1920

    - -

    -



    - -

    - - -

    Linear regression in Python

    -
    -
    -

    -

    - - -

    import numpy as np
-import matplotlib.pyplot as plt
-from IPython.display import display
-import sklearn
-from sklearn.linear_model import LinearRegression
-from sklearn.tree import DecisionTreeRegressor
-
-
-data = np.loadtxt('src/Hudson_Bay.csv', delimiter=',', skiprows=1)
-x = data[:,0]
-y = data[:,1]
-line = np.linspace(1900,1920,1000,endpoint=False).reshape(-1,1)
-reg = DecisionTreeRegressor(min_samples_split=3).fit(x.reshape(-1,1),y.reshape(-1,1))
-plt.plot(line, reg.predict(line), label="decision tree")
-regline = LinearRegression().fit(x.reshape(-1,1),y.reshape(-1,1))
-plt.plot(line, regline.predict(line), label= "Linear Regression")
-plt.plot(x, y, label= "Linear Regression")
-plt.show()
-
    -

    -

    -
    - - -

    - - -

    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. - -

    -This behaviour has inspired a simple mathematical model for an artificial neuron. - -$$ - y = f\left(\sum_{i=1}^n w_ix_i\right) = f(u) -$$ - -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. -

    -
    - - -

    - - -

    A simple perceptron model

    - -

    -



    - -

    - - -

    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. - -

    -



    -
    -
    - - -

    - - -

    The system: two electrons in a harmonic oscillator trap in two dimensions

    - -

    -The Hamiltonian of the quantum dot is given by -$$ \hat{H} = \hat{H}_0 + \hat{V}, -$$ - -where \( \hat{H}_0 \) is the many-body HO Hamiltonian, and \( \hat{V} \) is the -inter-electron Coulomb interactions. In dimensionless units, -$$ \hat{V}= \sum_{i < j}^N \frac{1}{r_{ij}}, -$$ - -with \( r_{ij}=\sqrt{\mathbf{r}_i^2 - \mathbf{r}_j^2} \). - -

    -This leads to the separable Hamiltonian, with the relative motion part given by (\( r_{ij}=r \)) -$$ -\hat{H}_r=-\nabla^2_r + \frac{1}{4}\omega^2r^2+ \frac{1}{r}, -$$ - -plus a standard Harmonic Oscillator problem for the center-of-mass motion. -This system has analytical solutions in two and three dimensions (M. Taut 1993 and 1994). - -

    - - -

    Quantum Monte Carlo Motivation

    -
    -
    -

    -Given a hamiltonian \( H \) and a trial wave function \( \Psi_T \), the variational principle states that the expectation value of \( \langle H \rangle \), defined through -$$ - \langle E \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})H(\boldsymbol{R})\Psi_T(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}, -$$ - -is an upper bound to the ground state energy \( E_0 \) of the hamiltonian \( H \), that is -$$ - E_0 \le \langle H \rangle . -$$ - -In general, the integrals involved in the calculation of various expectation values are multi-dimensional ones. Traditional integration methods such as the Gauss-Legendre will not be adequate for say the computation of the energy of a many-body system. -

    -
    - - -

    - - -

    Quantum Monte Carlo Motivation

    -
    -
    -

    -Choose a trial wave function -\( \psi_T(\boldsymbol{R}) \). -$$ - P(\boldsymbol{R},\boldsymbol{\alpha})= \frac{\left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2}{\int \left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2d\boldsymbol{R}}. -$$ - -This is our model, or likelihood/probability distribution function (PDF). It depends on some variational parameters \( \boldsymbol{\alpha} \). -The approximation to the expectation value of the Hamiltonian is now -$$ - \langle E[\boldsymbol{\alpha}] \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R},\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_T(\boldsymbol{R},\boldsymbol{\alpha})} - {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R},\boldsymbol{\alpha})\Psi_T(\boldsymbol{R},\boldsymbol{\alpha})}. -$$ -

    -
    - - -

    - - -

    Quantum Monte Carlo Motivation

    -
    -
    -

    -$$ - E_L(\boldsymbol{R},\boldsymbol{\alpha})=\frac{1}{\psi_T(\boldsymbol{R},\boldsymbol{\alpha})}H\psi_T(\boldsymbol{R},\boldsymbol{\alpha}), -$$ - -called the local energy, which, together with our trial PDF yields -$$ - E[\boldsymbol{\alpha}]=\int P(\boldsymbol{R})E_L(\boldsymbol{R},\boldsymbol{\alpha}) d\boldsymbol{R}\approx \frac{1}{N}\sum_{i=1}^NE_L(\boldsymbol{R_i},\boldsymbol{\alpha}) -$$ - -with \( N \) being the number of Monte Carlo samples. -

    -
    - - -

    - - -

    Quantum Monte Carlo

    -
    -
    -

    -The Algorithm for performing a variational Monte Carlo calculations runs thus as this - -

      -
    • Initialisation: Fix the number of Monte Carlo steps. Choose an initial \( \boldsymbol{R} \) and variational parameters \( \alpha \) and calculate \( \left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2 \).
      • -
    • Initialise the energy and the variance and start the Monte Carlo calculation by looping over trials.
      • - -
        -
      • Calculate a trial position \( \boldsymbol{R}_p=\boldsymbol{R}+r*step \) where \( r \) is a random variable \( r \in [0,1] \).
        • -
      • Metropolis algorithm to accept or reject this move \( w = P(\boldsymbol{R}_p,\boldsymbol{\alpha})/P(\boldsymbol{R},\boldsymbol{\alpha}) \).
        • -
      • If the step is accepted, then we set \( \boldsymbol{R}=\boldsymbol{R}_p \).
        • -
      • Update averages
        • -
      - -
    • Finish and compute final averages.
      • -
    - -Observe that the jumping in space is governed by the variable step. This is often called brute-force sampling. -Need importance sampling to get more relevant sampling. -
    -
    - - -

    - - -

    The trial wave function

    -
    -
    -

    -We want to perform a Variational Monte Carlo calculation of the ground state of two electrons in a quantum dot well with different oscillator energies, assuming total spin \( S=0 \). -Our trial wave function has the following form -$$ -\begin{equation} - \psi_{T}(\boldsymbol{r}_1,\boldsymbol{r}_2) = - C\exp{\left(-\alpha_1\omega(r_1^2+r_2^2)/2\right)} - \exp{\left(\frac{r_{12}}{(1+\alpha_2 r_{12})}\right)}, -\label{eq:trial} -\end{equation} -$$ - -where the $\alpha$s represent our variational parameters, two in this case. - -

    -Why does the trial function look like this? How did we get there? This will be our main motivation for switching to -Machine Learning. - -

    -

    -
    - - -

    - - -

    The correlation part of the wave function

    - -

    -To find an ansatz for the correlated part of the wave function, it is useful to rewrite the two-particle -local energy in terms of the relative and center-of-mass motion. -Let us denote the distance between the two electrons as -\( r_{12} \). We omit the center-of-mass motion since we are only interested in the case when -\( r_{12} \rightarrow 0 \). The contribution from the center-of-mass (CoM) variable \( {\bf R}_{\mathrm{CoM}} \) -gives only a finite contribution. -We focus only on the terms that are relevant for \( r_{12} \) and for three dimensions. The relevant local energy becomes then -$$ -\lim_{r_{12} \rightarrow 0}E_L(R)= - \frac{1}{{\cal R}_T(r_{12})}\left(2\frac{d^2}{dr_{ij}^2}+\frac{4}{r_{ij}}\frac{d}{dr_{ij}}+ -\frac{2}{r_{ij}}-\frac{l(l+1)}{r_{ij}^2}+2E -\right){\cal R}_T(r_{12}) = 0. -$$ - -Set \( l=0 \) and we have the so-called cusp condition -$$ -\frac{d {\cal R}_T(r_{12})}{dr_{12}} = -\frac{1}{2(l+1)} -{\cal R}_T(r_{12})\qquad r_{12}\to 0 -$$ - -

    - - -

    Resulting ansatz

    -The above results in -$$ -{\cal R}_T \propto \exp{(r_{ij}/2)}, -$$ - -for anti-parallel spins and -$$ -{\cal R}_T \propto \exp{(r_{ij}/4)}, -$$ - -for anti-parallel spins. -This is the so-called cusp condition for the relative motion, resulting in a minimal requirement -for the correlation part of the wave fuction. -For general systems containing more than say two electrons, we have this -condition for each electron pair \( ij \). - -

    - - -

    The VMC code

    - -

    - - -

    # Importing various packages
-from math import exp, sqrt
-from random import random, seed
-import numpy as np
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-from matplotlib import cm
-from matplotlib.ticker import LinearLocator, FormatStrFormatter
-import sys
-
-#Trial wave function for quantum dots in two dims
-def WaveFunction(r,alpha,beta):
-    r1 = r[0,0]**2 + r[0,1]**2
-    r2 = r[1,0]**2 + r[1,1]**2
-    r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
-    deno = r12/(1+beta*r12)
-    return exp(-0.5*alpha*(r1+r2)+deno)
-
-#Local energy  for quantum dots in two dims, using analytical local energy
-def LocalEnergy(r,alpha,beta):
-    
-    r1 = (r[0,0]**2 + r[0,1]**2)
-    r2 = (r[1,0]**2 + r[1,1]**2)
-    r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
-    deno = 1.0/(1+beta*r12)
-    deno2 = deno*deno
-    return 0.5*(1-alpha*alpha)*(r1 + r2) +2.0*alpha + 1.0/r12+deno2*(alpha*r12-deno2+2*beta*deno-1.0/r12)
-
-# The Monte Carlo sampling with the Metropolis algo
-def MonteCarloSampling():
-
-    NumberMCcycles= 100000
-    StepSize = 1.0
-    # positions
-    PositionOld = np.zeros((NumberParticles,Dimension), np.double)
-    PositionNew = np.zeros((NumberParticles,Dimension), np.double)
-    # seed for rng generator
-    seed()
-    # start variational parameter
-    alpha = 0.9
-    for ia in range(MaxVariations):
-        alpha += .025
-        AlphaValues[ia] = alpha
-        beta = 0.2 
-        for jb in range(MaxVariations):
-            beta += .01
-            BetaValues[jb] = beta
-            energy = energy2 = 0.0
-            DeltaE = 0.0
-            #Initial position
-            for i in range(NumberParticles):
-                for j in range(Dimension):
-                    PositionOld[i,j] = StepSize * (random() - .5)
-            wfold = WaveFunction(PositionOld,alpha,beta)
-
-            #Loop over MC MCcycles
-            for MCcycle in range(NumberMCcycles):
-                #Trial position
-                for i in range(NumberParticles):
-                    for j in range(Dimension):
-                        PositionNew[i,j] = PositionOld[i,j] + StepSize * (random() - .5)
-                wfnew = WaveFunction(PositionNew,alpha,beta)
-
-                #Metropolis test to see whether we accept the move
-                if random() < wfnew**2 / wfold**2:
-                   PositionOld = PositionNew.copy()
-                   wfold = wfnew
-                   DeltaE = LocalEnergy(PositionOld,alpha,beta)
-                energy += DeltaE
-                energy2 += DeltaE**2
-
-            #We calculate mean, variance and error ...
-            energy /= NumberMCcycles
-            energy2 /= NumberMCcycles
-            variance = energy2 - energy**2
-            error = sqrt(variance/NumberMCcycles)
-            Energies[ia,jb] = energy    
-    return Energies, AlphaValues, BetaValues
-
-
-#Here starts the main program with variable declarations
-NumberParticles = 2
-Dimension = 2
-MaxVariations = 10
-Energies = np.zeros((MaxVariations,MaxVariations))
-AlphaValues = np.zeros(MaxVariations)
-BetaValues = np.zeros(MaxVariations)
-(Energies, AlphaValues, BetaValues) = MonteCarloSampling()
-
-# Prepare for plots
-fig = plt.figure()
-ax = fig.gca(projection='3d')
-# Plot the surface.
-X, Y = np.meshgrid(AlphaValues, BetaValues)
-surf = ax.plot_surface(X, Y, Energies,cmap=cm.coolwarm,linewidth=0, antialiased=False)
-# Customize the z axis.
-zmin = np.matrix(Energies).min()
-zmax = np.matrix(Energies).max()
-ax.set_zlim(zmin, zmax)
-ax.set_xlabel(r'$\alpha$')
-ax.set_ylabel(r'$\beta$')
-ax.set_zlabel(r'$\langle E \rangle$')
-ax.zaxis.set_major_locator(LinearLocator(10))
-ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
-# Add a color bar which maps values to colors.
-fig.colorbar(surf, shrink=0.5, aspect=5)
-plt.show()
-
    -

    - - -

    Technical aspect, improvements and how to define the cost function

    -
    -
    -

    - -

    -The above procedure is not the smartest one. Looping over all variational parameters becomes expensive. -Also, we don't use importance sampling and optimizations of the standard deviation (blocking, bootstrap, jackknife). -Such codes are included in the above Github address. - -

    -We can also be smarter and use minimization methods to find the optimal variational parameters with fewer Monte Carlo cycles and then -fire up our heavy artillery. - -

    -One way to achieve this is to minimize the energy as function of the variational parameters. - -

    -

    -
    - - -

    - - -

    Energy derivatives

    -
    -
    -

    -To find the derivatives of the local energy expectation value as function of the variational parameters, we can use the chain rule and the hermiticity of the Hamiltonian. - -

    -Let us define -$$ -\bar{E}_{\alpha_i}=\frac{d\langle E_L\rangle}{d\alpha_i}. -$$ - -as the derivative of the energy with respect to the variational parameter \( \alpha_i \) -We define also the derivative of the trial function (skipping the subindex \( T \)) as -$$ -\bar{\Psi}_{i}=\frac{d\Psi}{d\alpha_i}. -$$ -

    -
    - - -

    - - -

    Derivatives of the local energy

    -
    -
    -

    -The elements of the gradient of the local energy are then (using the chain rule and the hermiticity of the Hamiltonian) -$$ -\bar{E}_{i}= 2\left( \langle \frac{\bar{\Psi}_{i}}{\Psi}E_L\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L \rangle\right). -$$ - -From a computational point of view it means that you need to compute the expectation values of -$$ -\langle \frac{\bar{\Psi}_{i}}{\Psi}E_L\rangle, -$$ - -and -$$ -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L\rangle -$$ - -These integrals are evaluted using MC intergration (with all its possible error sources). -We can then use methods like stochastic gradient or other minimization methods to find the optimal variational parameters (I don't discuss this topic here, but these methods are very important in ML). -

    -
    - - -

    - - -

    How do we define our cost function?

    -
    -
    -

    -We have a model, our likelihood function. - -

    -How should we define the cost function? - -

    -

    -
    - - -

    - - -

    Meet the variance and its derivatives

    -
    -
    -

    -Suppose the trial function (our model) is the exact wave function. The action of the hamiltionan -on the wave function -$$ - H\Psi = \mathrm{constant}\times \Psi, -$$ - -The integral which defines various -expectation values involving moments of the hamiltonian becomes then -$$ - \langle E^n \rangle = \langle H^n \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})H^n(\boldsymbol{R})\Psi(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})}= -\mathrm{constant}\times\frac{\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})}=\mathrm{constant}. -$$ - -This gives an important information: If I want the variance, the exact wave function leads to zero variance! -The variance is defined as -$$ -\sigma_E = \langle E^2\rangle - \langle E\rangle^2. -$$ - -Variation is then performed by minimizing both the energy and the variance. - -

    -

    -
    - - -

    - - -

    The variance defines the cost function

    -
    -
    -

    -We can then take the derivatives of -$$ -\sigma_E = \langle E^2\rangle - \langle E\rangle^2, -$$ - -with respect to the variational parameters. The derivatives of the variance can then be used to defined the -so-called Hessian matrix, which in turn allows us to use minimization methods like Newton's method or -standard gradient methods. - -

    -This leads to however a more complicated expression, with obvious errors when evaluating integrals by Monte Carlo integration. Less used, see however Filippi and Umrigar. The expression becomes complicated -$$ -\bar{E}_{ij} = 2\left[ \langle (\frac{\bar{\Psi}_{ij}}{\Psi}+\frac{\bar{\Psi}_{j}}{\Psi}\frac{\bar{\Psi}_{i}}{\Psi})(E_L-\langle E\rangle)\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\bar{E}_j-\langle \frac{\bar{\Psi}_{j}}{\Psi}\rangle\bar{E}_i\right] +\langle \frac{\bar{\Psi}_{i}}{\Psi}E_L{_j}\rangle +\langle \frac{\bar{\Psi}_{j}}{\Psi}E_L{_i}\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L{_j}\rangle \langle \frac{\bar{\Psi}_{j}}{\Psi}\rangle\langle E_L{_i}\rangle. -$$ - -

    -Evaluating the cost function means having to evaluate the above second derivative of the energy. - -

    -

    -
    - - -

    - - -

    Boltzmann Machines

    +

    Boltzmann Machines

    Why use a generative model rather than the more well known discriminative deep neural networks (DNN)? @@ -1146,7 +243,7 @@ Why use a generative model rather than the more well known discriminative deep n -

    Some similarities and differences from DNNs

    +

    Some similarities and differences from DNNs

    1. Both use gradient-descent based learning procedures for minimizing cost functions
      2. @@ -1159,7 +256,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some

      -

      The structure of the RBM network

      +

      The structure of the RBM network





      @@ -1167,7 +264,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some

      -

      The network

      +

      The network

      The network layers: @@ -1179,7 +276,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some -

      Goals

      +

      Goals

      The goal of the hidden layer is to increase the model's expressive power. We 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). @@ -1203,7 +300,7 @@ Examples of this trick being employed in physics: -

      Joint distribution and the Energy function

      +

      Joint distribution and the Energy function

      The restricted Boltzmann machine is described by a Bolztmann distribution $$ \begin{align} @@ -1225,7 +322,7 @@ It is common to ignore \( T_0 \) by setting it to one.

      -

      Network Elements

      +

      Network Elements

      The function \( E(\mathbf{x},\mathbf{h}) \) gives the energy of a @@ -1238,7 +335,7 @@ adjusting the energy function to best fit our problem.

      -

      Defining different types of RBMs

      +

      Defining different types of RBMs

      There are different variants of RBMs, and the differences lie in the types of visible and hidden units we choose as well as in the implementation of the energy function \( E(\mathbf{x},\mathbf{h}) \).

      @@ -1278,7 +375,7 @@ $$

      -

      More about RBMs

      +

      More about RBMs

      1. Useful when we model continuous data (i.e., we wish \( \mathbf{x} \) to be continuous)
        2. @@ -1296,7 +393,7 @@ Other types of units include: -

        Sampling: Metropolis sampling

        +

        Sampling: Metropolis sampling

        In order to sample from the RBM probability distribution it is common to use Markov Chain Monte Carlo (MCMC) algorithms such as Metropolis-Hastings or Gibbs sampling.

        @@ -1321,7 +418,7 @@ $$

        -

        Sampling: Gibbs sampling

        +

        Sampling: Gibbs sampling

        In this method we sample from the joint probability \( P_{rbm} (\mathbf{x}, \mathbf{h}) \) by way of a two step sampling process. We alternately update the visible and hidden units. @@ -1343,7 +440,7 @@ where we recognize the logistic sigmoid function \( \sigma (x) = 1/(1+exp(-x)) \

        -

        Gaussian RBM

        +

        Gaussian RBM

        For the Gaussian-Binary RBM the conditional probabilities are $$ \begin{align} @@ -1359,7 +456,7 @@ while the visible units now follow a normal distribution, we see the hidden unit

        -

        Cost function

        +

        Cost function

        When working with a training dataset, the most common training approach is maximizing the log-likelihood of the training data. The log likelihood characterizes the log-probability of generating the observed data using our generative model. Using this method our cost function is chosen as the negative log-likelihood. The learning then consists of trying to find parameters that maximize the probability of the dataset, and is known as Maximum Likelihood Estimation (MLE). @@ -1379,7 +476,7 @@ Our cost function is the negative log-likelihood, \( \mathcal{C}(\{ \theta_i \})

        -

        Optimization / Training

        +

        Optimization / Training

        The training procedure of choice often is Stochastic Gradient Descent (SGD). It consists of a series of iterations where we update the parameters according to the equation $$ \begin{align} @@ -1422,7 +519,7 @@ $$

        -

        More on RBMs

        +

        More on RBMs

        The data-dependent term in the gradient is known as the positive phase of the gradient, while the model-dependent term is known as the negative phase of the gradient. The aim of the training is to lower the energy of configurations that are near observed data points (increasing their probability), and raising the energy of configurations that are far from observed data points (decreasing their probability). @@ -1446,7 +543,7 @@ To get the expecation values with respect to the data, we set the visib

        -

        Which sampling to use

        +

        Which sampling to use

        To get the expectation values with respect to the model, we use Gibbs sampling. We can either initialize the \( \boldsymbol{x} \) randomly or with a training sample. While we ideally want a large number of Gibbs iterations \( n\rightarrow n \), one might decide to truncate it earlier for efficiency. Doing this while having intialized \( \boldsymbol{x} \) with a training data vector is referred to as contrastive divergence (CD), because one is then closer to approximating the gradient of this function than the negative log-likelihood. The contrastive divergence function is the difference between two Kullback-Leibler divergences (also called relative entropy), which measure how one probability distribution diverges from a second, expected probability distribution (in this case the estimated one from the ground truth one). @@ -1454,7 +551,7 @@ To get the expectation values with respect to the model, we use Gibbs s

        -

        RBMs for the quantum many body problem

        +

        RBMs for the quantum many body problem

        The idea of applying RBMs to quantum many body problems was presented by G. Carleo and M. Troyer, working with ETH Zurich and Microsoft Research. @@ -1478,7 +575,7 @@ Some of their motivation included -

        Choose the right RBM

        +

        Choose the right RBM

        Carleo and Troyer applied the RBM to the quantum mechanical spin lattice systems of the Ising model and Heisenberg model, with encouraging results. Our goal is to test the method on systems of moving particles. For the spin lattice systems it was natural to use a binary-binary RBM, with the nodes taking values of 1 and -1. For moving particles, on the other hand, we want the visible nodes to be continuous, representing position coordinates. Thus, we start by choosing a Gaussian-binary RBM, where the visible nodes are continuous and hidden nodes take on values of 0 or 1. If eventually we would like the hidden nodes to be continuous as well the rectified linear units seem like the most relevant choice. @@ -1486,7 +583,7 @@ Carleo and Troyer applied the RBM to the quantum mechanical spin lattice systems

        -

        Representing the wave function

        +

        Representing the wave function

        The wavefunction should be a probability amplitude depending on \( \boldsymbol{x} \). The RBM model is given by the joint distribution of \( \boldsymbol{x} \) and \( \boldsymbol{h} \) $$ \begin{align} @@ -1523,7 +620,7 @@ $$

        -

        Choose the cost function

        +

        Choose the cost function

        Now we don't necessarily have training data (unless we generate it by using some other method). However, what we do have is the variational principle which allows us to obtain the ground state wave function by minimizing the expectation value of the energy of a trial wavefunction (corresponding to the untrained NQS). Similarly to the traditional variational Monte Carlo method then, it is the local energy we wish to minimize. The gradient to use for the stochastic gradient descent procedure is $$ \begin{align} @@ -1544,7 +641,7 @@ $$

        -

        Running the codes

        +

        Running the codes

        @@ -1564,7 +661,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat

        -

        Energy as function of iterations, \( N=2 \) electrons

        +

        Energy as function of iterations, \( N=2 \) electrons

        @@ -1576,7 +673,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat

        -

        Energy as function of iterations, \( N=6 \) electrons

        +

        Energy as function of iterations, \( N=6 \) electrons

        @@ -1588,7 +685,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat

        -

        Conclusions and where do we stand

        +

        Conclusions and where do we stand

        diff --git a/doc/pub/BM/html/BM-reveal.html b/doc/pub/BM/html/BM-reveal.html index d84d03de9..66b338087 100644 --- a/doc/pub/BM/html/BM-reveal.html +++ b/doc/pub/BM/html/BM-reveal.html @@ -3,9 +3,9 @@ - + -Solving quantum mechanical problems with Machine Learning and Boltzmann machines +Machine Learning and Boltzmann machines @@ -132,7 +132,7 @@ MathJax.Hub.Config({ -

        Solving quantum mechanical problems with Machine Learning and Boltzmann machines

        +

        Machine Learning and Boltzmann machines

        @@ -147,7 +147,7 @@ MathJax.Hub.Config({

        Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University and Department of Physics, University of Oslo, Norway

         
        -

        Jul 26, 2018

        +

        Nov 19, 2018


        @@ -158,198 +158,16 @@ MathJax.Hub.Config({

        -

        What is this talk about?

        +

        Unsupervised learning, ovrarching aims

        -The main aim is to give you a short and pedestrian introduction to how we can use Machine Learning methods -to solve quantum mechanical many-body problems. And why this could be of interest. - -

        -The hope is that after this talk you have gotten the basic ideas to get you started. Peeping into https://github.com/mhjensenseminars/MachineLearningTalk, you'll find a Jupyter notebook, slides, codes etc that will allow you to reproduce the simulations discussed here, and perhaps run your own very first calculations. - -

        -Furthermore, since we will also talk about Monte Carlo simulations, the exciting idea of Monte Carlo optimization and Machine Learning can be used to mesh strategies for reducing the curse of dimensionality. - -

        -

        Why?

        - -

        -How can we avoid the dimensionality curse? Many possibilities - -

          -

        1. smarter basis functions
          2. -

        2. resummation of specific correlations
          3. -

        3. stochastic sampling of high-lying states (stochastic FCI, CC and soon SRG)
          4. -

        4. many more
          5. -
        -

        - -Machine Learning and Quantum Computing hold also great promise in tackling the -ever increasing dimensionalities. Here we will focus on Machine Learning. -

        - - -
        -

        Overview

        -
        - -
          -

        • Short intro to Machine Learning
          • -

        • Variational Monte Carlo (Markov Chain Monte Carlo, \( \mathrm{MC}^2 \)) and two-electron quantum dots, solving quantum mechanical problems in a stochastic way. It will serve as our motivation for switching to Machine Learning.
          • -

        • From Variational Monte Carlo to Boltzmann Machines and Machine Learning
          • -
        -
        -
        - - -
        -

        What are the Machine Learning calculations here based on?

        - -

        -This work is inspired by the idea of representing the wave function with -a restricted Boltzmann machine (RBM), presented recently by G. Carleo and M. Troyer, Science 355, Issue 6325, pp. 602-606 (2017). They -named such a wave function/network a neural network quantum state (NQS). In their article they apply it to the quantum mechanical -spin lattice systems of the Ising model and Heisenberg model, with -encouraging results. - -

        -Thanks to Vilde Flugsrud and Alfred Alocias Mariadason for many discussions and interpretations of the above paper. -

        - - -
        -

        A new world

        -
        - -

        -Machine learning (ML) 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 ML are - -

        -

        - -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. -

        -
        - - -
        -

        Lots of room for creativity

        -
        - -

        -Not all the -algorithms and methods can be given a rigorous mathematical -justification, opening up thereby for experimenting -and trial and error and thereby exciting new developments. -

        - -

        -

        - -

        -A solid command of linear algebra, multivariate theory, -probability theory, statistical data analysis, -understanding errors and Monte Carlo methods is important in order to understand many of the -various algorithms and methods. -

        - -

        -

        -Job market, a personal statement. -

        -A familiarity with ML is almost becoming a prerequisite for many of the most exciting employment -opportunities. Add quantum computing and there you are. -

        -
        - - -
        -

        Some members of the ML family

        - -

        -

        - -
          -

        1. Linear regression and its variants, in essence polynomial regression
          2. -

        2. Decision tree algorithms, from simpler to more complex ones like random forests
          3. -

        3. Nearest neighbors models
          4. -

        4. Bayesian statistics and regression
          5. -

        5. Support vector machines and finally various variants of
          6. -

        6. Artifical neural networks and deep learning
          7. -

        7. ....
          8. -
        -
        -
        - - -
        -

        What are the basic ingredients?

        -
        - -

        -Almost every problem in ML and data science starts with the same ingredients: - -

          -

        • The dataset \( \mathbf{x} \) (could be some observable quantity of the system we are studying)
          • -

        • A model which is a function of a set of parameters \( \mathbf{\alpha} \) that relates to the dataset, say a likelihood function \( p(\mathbf{x}\vert \mathbf{\alpha}) \) or just a simple model \( f(\mathbf{\alpha}) \)
          • -

        • A so-called cost function \( \mathcal{C} (\mathbf{x}, f(\mathbf{\alpha})) \) which allows us to decide how well our model represents the dataset.
          • -
        -

        - -We seek to minimize the function \( \mathcal{C} (\mathbf{x}, f(\mathbf{\alpha})) \) by finding the parameter values which minimize \( \mathcal{C} \). This leads to various minimization algorithms. -

        -
        - - -
        -

        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 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. - -

        -You will notice however that many of the basic ideas discussed do come from Physics! -

        -
        - - -
        -

        Types of Machine Learning

        +

        Types of Machine Learning

        @@ -380,648 +198,7 @@ Some of the most common tasks are:
        -

        References

        -
        - - -
        -
        - - -
        -

        Another interesting article

        - -

        -Here we will use so-called reduced Boltzmann Machines to simulate quantum many-body problems. For Monte Carlo aficionados, there is a very close similarity with what are called shadow wave functions, see the work of Pederiva and Kalos and collaborators, Phys Rev. E 90, 053304 (2014). -

        - - -
        -

        Just for fun illustration of scikit-learn: Hudson bay data

        -
        - -

        -Most mammalian predators rely on a variety of prey, which complicates mathematical modeling; however, a few predators have become highly specialized and seek almost exclusively a single prey species. An example of this simplified predator-prey interaction is seen in Canadian northern forests, where the populations of the lynx and the snowshoe hare are intertwined in a life and death struggle. - -

        -One reason that this particular system has been so extensively studied is that the Hudson Bay company kept careful records of all furs from the early 1800s into the 1900s. The records for the furs collected by the Hudson Bay company showed distinct oscillations (approximately 12 year periods), suggesting that these species caused almost periodic fluctuations of each other's populations. The table here shows data from 1900 to 1920. - - -

        -
        - - -
        -

        Hares and lynx in Hudson bay from 1900 to 1920

        - -

        -



        -
        - - -
        -

        Linear regression in Python

        -
        - -

        - - -

        import numpy as np
-import matplotlib.pyplot as plt
-from IPython.display import display
-import sklearn
-from sklearn.linear_model import LinearRegression
-from sklearn.tree import DecisionTreeRegressor
-
-
-data = np.loadtxt('src/Hudson_Bay.csv', delimiter=',', skiprows=1)
-x = data[:,0]
-y = data[:,1]
-line = np.linspace(1900,1920,1000,endpoint=False).reshape(-1,1)
-reg = DecisionTreeRegressor(min_samples_split=3).fit(x.reshape(-1,1),y.reshape(-1,1))
-plt.plot(line, reg.predict(line), label="decision tree")
-regline = LinearRegression().fit(x.reshape(-1,1),y.reshape(-1,1))
-plt.plot(line, regline.predict(line), label= "Linear Regression")
-plt.plot(x, y, label= "Linear Regression")
-plt.show()
-
        - -
        -
        - - -
        -

        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. - -

        -This behaviour has inspired a simple mathematical model for an artificial neuron. - -

         
        -$$ - y = f\left(\sum_{i=1}^n w_ix_i\right) = f(u) -$$ -

         
        - -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. -

        -
        - - -
        -

        A simple perceptron model

        - -

        -



        -
        - - -
        -

        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. - -

        -



        -
        -
        - - -
        -

        The system: two electrons in a harmonic oscillator trap in two dimensions

        - -

        -The Hamiltonian of the quantum dot is given by -

         
        -$$ \hat{H} = \hat{H}_0 + \hat{V}, -$$ -

         
        - -where \( \hat{H}_0 \) is the many-body HO Hamiltonian, and \( \hat{V} \) is the -inter-electron Coulomb interactions. In dimensionless units, -

         
        -$$ \hat{V}= \sum_{i < j}^N \frac{1}{r_{ij}}, -$$ -

         
        - -with \( r_{ij}=\sqrt{\mathbf{r}_i^2 - \mathbf{r}_j^2} \). - -

        -This leads to the separable Hamiltonian, with the relative motion part given by (\( r_{ij}=r \)) -

         
        -$$ -\hat{H}_r=-\nabla^2_r + \frac{1}{4}\omega^2r^2+ \frac{1}{r}, -$$ -

         
        - -plus a standard Harmonic Oscillator problem for the center-of-mass motion. -This system has analytical solutions in two and three dimensions (M. Taut 1993 and 1994). -

        - - -
        -

        Quantum Monte Carlo Motivation

        -
        - -

        -Given a hamiltonian \( H \) and a trial wave function \( \Psi_T \), the variational principle states that the expectation value of \( \langle H \rangle \), defined through -

         
        -$$ - \langle E \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})H(\boldsymbol{R})\Psi_T(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}, -$$ -

         
        - -is an upper bound to the ground state energy \( E_0 \) of the hamiltonian \( H \), that is -

         
        -$$ - E_0 \le \langle H \rangle . -$$ -

         
        - -In general, the integrals involved in the calculation of various expectation values are multi-dimensional ones. Traditional integration methods such as the Gauss-Legendre will not be adequate for say the computation of the energy of a many-body system. -

        -
        - - -
        -

        Quantum Monte Carlo Motivation

        -
        -Basic steps. -

        -Choose a trial wave function -\( \psi_T(\boldsymbol{R}) \). -

         
        -$$ - P(\boldsymbol{R},\boldsymbol{\alpha})= \frac{\left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2}{\int \left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2d\boldsymbol{R}}. -$$ -

         
        - -This is our model, or likelihood/probability distribution function (PDF). It depends on some variational parameters \( \boldsymbol{\alpha} \). -The approximation to the expectation value of the Hamiltonian is now -

         
        -$$ - \langle E[\boldsymbol{\alpha}] \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R},\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_T(\boldsymbol{R},\boldsymbol{\alpha})} - {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R},\boldsymbol{\alpha})\Psi_T(\boldsymbol{R},\boldsymbol{\alpha})}. -$$ -

         
        -

        -
        - - -
        -

        Quantum Monte Carlo Motivation

        -
        -Define a new quantity. -

         
        -$$ - E_L(\boldsymbol{R},\boldsymbol{\alpha})=\frac{1}{\psi_T(\boldsymbol{R},\boldsymbol{\alpha})}H\psi_T(\boldsymbol{R},\boldsymbol{\alpha}), -$$ -

         
        - -called the local energy, which, together with our trial PDF yields -

         
        -$$ - E[\boldsymbol{\alpha}]=\int P(\boldsymbol{R})E_L(\boldsymbol{R},\boldsymbol{\alpha}) d\boldsymbol{R}\approx \frac{1}{N}\sum_{i=1}^NE_L(\boldsymbol{R_i},\boldsymbol{\alpha}) -$$ -

         
        - -with \( N \) being the number of Monte Carlo samples. -

        -
        - - -
        -

        Quantum Monte Carlo

        -
        - -

        -The Algorithm for performing a variational Monte Carlo calculations runs thus as this - -

          - -

        • Initialisation: Fix the number of Monte Carlo steps. Choose an initial \( \boldsymbol{R} \) and variational parameters \( \alpha \) and calculate \( \left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2 \).
          • - -

        • Initialise the energy and the variance and start the Monte Carlo calculation by looping over trials.
          • - -
            - -

          • Calculate a trial position \( \boldsymbol{R}_p=\boldsymbol{R}+r*step \) where \( r \) is a random variable \( r \in [0,1] \).
            • - -

          • Metropolis algorithm to accept or reject this move \( w = P(\boldsymbol{R}_p,\boldsymbol{\alpha})/P(\boldsymbol{R},\boldsymbol{\alpha}) \).
            • - -

          • If the step is accepted, then we set \( \boldsymbol{R}=\boldsymbol{R}_p \).
            • - -

          • Update averages
            • -
          -

        • Finish and compute final averages.
          • -
        -

        - -Observe that the jumping in space is governed by the variable step. This is often called brute-force sampling. -Need importance sampling to get more relevant sampling. -

        -
        - - -
        -

        The trial wave function

        -
        - -

        -We want to perform a Variational Monte Carlo calculation of the ground state of two electrons in a quantum dot well with different oscillator energies, assuming total spin \( S=0 \). -Our trial wave function has the following form -

         
        -$$ -\begin{equation} - \psi_{T}(\boldsymbol{r}_1,\boldsymbol{r}_2) = - C\exp{\left(-\alpha_1\omega(r_1^2+r_2^2)/2\right)} - \exp{\left(\frac{r_{12}}{(1+\alpha_2 r_{12})}\right)}, -\tag{1} -\end{equation} -$$ -

         
        - -where the $\alpha$s represent our variational parameters, two in this case. - -

        -Why does the trial function look like this? How did we get there? This will be our main motivation for switching to -Machine Learning. - - -

        -
        - - -
        -

        The correlation part of the wave function

        - -

        -To find an ansatz for the correlated part of the wave function, it is useful to rewrite the two-particle -local energy in terms of the relative and center-of-mass motion. -Let us denote the distance between the two electrons as -\( r_{12} \). We omit the center-of-mass motion since we are only interested in the case when -\( r_{12} \rightarrow 0 \). The contribution from the center-of-mass (CoM) variable \( {\bf R}_{\mathrm{CoM}} \) -gives only a finite contribution. -We focus only on the terms that are relevant for \( r_{12} \) and for three dimensions. The relevant local energy becomes then -

         
        -$$ -\lim_{r_{12} \rightarrow 0}E_L(R)= - \frac{1}{{\cal R}_T(r_{12})}\left(2\frac{d^2}{dr_{ij}^2}+\frac{4}{r_{ij}}\frac{d}{dr_{ij}}+ -\frac{2}{r_{ij}}-\frac{l(l+1)}{r_{ij}^2}+2E -\right){\cal R}_T(r_{12}) = 0. -$$ -

         
        - -Set \( l=0 \) and we have the so-called cusp condition -

         
        -$$ -\frac{d {\cal R}_T(r_{12})}{dr_{12}} = -\frac{1}{2(l+1)} -{\cal R}_T(r_{12})\qquad r_{12}\to 0 -$$ -

         
        -

        - - -
        -

        Resulting ansatz

        -The above results in -

         
        -$$ -{\cal R}_T \propto \exp{(r_{ij}/2)}, -$$ -

         
        - -for anti-parallel spins and -

         
        -$$ -{\cal R}_T \propto \exp{(r_{ij}/4)}, -$$ -

         
        - -for anti-parallel spins. -This is the so-called cusp condition for the relative motion, resulting in a minimal requirement -for the correlation part of the wave fuction. -For general systems containing more than say two electrons, we have this -condition for each electron pair \( ij \). -

        - - -
        -

        The VMC code

        - -

        - - -

        # Importing various packages
-from math import exp, sqrt
-from random import random, seed
-import numpy as np
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-from matplotlib import cm
-from matplotlib.ticker import LinearLocator, FormatStrFormatter
-import sys
-
-#Trial wave function for quantum dots in two dims
-def WaveFunction(r,alpha,beta):
-    r1 = r[0,0]**2 + r[0,1]**2
-    r2 = r[1,0]**2 + r[1,1]**2
-    r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
-    deno = r12/(1+beta*r12)
-    return exp(-0.5*alpha*(r1+r2)+deno)
-
-#Local energy  for quantum dots in two dims, using analytical local energy
-def LocalEnergy(r,alpha,beta):
-    
-    r1 = (r[0,0]**2 + r[0,1]**2)
-    r2 = (r[1,0]**2 + r[1,1]**2)
-    r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
-    deno = 1.0/(1+beta*r12)
-    deno2 = deno*deno
-    return 0.5*(1-alpha*alpha)*(r1 + r2) +2.0*alpha + 1.0/r12+deno2*(alpha*r12-deno2+2*beta*deno-1.0/r12)
-
-# The Monte Carlo sampling with the Metropolis algo
-def MonteCarloSampling():
-
-    NumberMCcycles= 100000
-    StepSize = 1.0
-    # positions
-    PositionOld = np.zeros((NumberParticles,Dimension), np.double)
-    PositionNew = np.zeros((NumberParticles,Dimension), np.double)
-    # seed for rng generator
-    seed()
-    # start variational parameter
-    alpha = 0.9
-    for ia in range(MaxVariations):
-        alpha += .025
-        AlphaValues[ia] = alpha
-        beta = 0.2 
-        for jb in range(MaxVariations):
-            beta += .01
-            BetaValues[jb] = beta
-            energy = energy2 = 0.0
-            DeltaE = 0.0
-            #Initial position
-            for i in range(NumberParticles):
-                for j in range(Dimension):
-                    PositionOld[i,j] = StepSize * (random() - .5)
-            wfold = WaveFunction(PositionOld,alpha,beta)
-
-            #Loop over MC MCcycles
-            for MCcycle in range(NumberMCcycles):
-                #Trial position
-                for i in range(NumberParticles):
-                    for j in range(Dimension):
-                        PositionNew[i,j] = PositionOld[i,j] + StepSize * (random() - .5)
-                wfnew = WaveFunction(PositionNew,alpha,beta)
-
-                #Metropolis test to see whether we accept the move
-                if random() < wfnew**2 / wfold**2:
-                   PositionOld = PositionNew.copy()
-                   wfold = wfnew
-                   DeltaE = LocalEnergy(PositionOld,alpha,beta)
-                energy += DeltaE
-                energy2 += DeltaE**2
-
-            #We calculate mean, variance and error ...
-            energy /= NumberMCcycles
-            energy2 /= NumberMCcycles
-            variance = energy2 - energy**2
-            error = sqrt(variance/NumberMCcycles)
-            Energies[ia,jb] = energy    
-    return Energies, AlphaValues, BetaValues
-
-
-#Here starts the main program with variable declarations
-NumberParticles = 2
-Dimension = 2
-MaxVariations = 10
-Energies = np.zeros((MaxVariations,MaxVariations))
-AlphaValues = np.zeros(MaxVariations)
-BetaValues = np.zeros(MaxVariations)
-(Energies, AlphaValues, BetaValues) = MonteCarloSampling()
-
-# Prepare for plots
-fig = plt.figure()
-ax = fig.gca(projection='3d')
-# Plot the surface.
-X, Y = np.meshgrid(AlphaValues, BetaValues)
-surf = ax.plot_surface(X, Y, Energies,cmap=cm.coolwarm,linewidth=0, antialiased=False)
-# Customize the z axis.
-zmin = np.matrix(Energies).min()
-zmax = np.matrix(Energies).max()
-ax.set_zlim(zmin, zmax)
-ax.set_xlabel(r'$\alpha$')
-ax.set_ylabel(r'$\beta$')
-ax.set_zlabel(r'$\langle E \rangle$')
-ax.zaxis.set_major_locator(LinearLocator(10))
-ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
-# Add a color bar which maps values to colors.
-fig.colorbar(surf, shrink=0.5, aspect=5)
-plt.show()
-
        -
        - - -
        -

        Technical aspect, improvements and how to define the cost function

        -
        - -

        -The above procedure is not the smartest one. Looping over all variational parameters becomes expensive. -Also, we don't use importance sampling and optimizations of the standard deviation (blocking, bootstrap, jackknife). -Such codes are included in the above Github address. - -

        -We can also be smarter and use minimization methods to find the optimal variational parameters with fewer Monte Carlo cycles and then -fire up our heavy artillery. - -

        -One way to achieve this is to minimize the energy as function of the variational parameters. - - -

        -
        - - -
        -

        Energy derivatives

        -
        - -

        -To find the derivatives of the local energy expectation value as function of the variational parameters, we can use the chain rule and the hermiticity of the Hamiltonian. - -

        -Let us define -

         
        -$$ -\bar{E}_{\alpha_i}=\frac{d\langle E_L\rangle}{d\alpha_i}. -$$ -

         
        - -as the derivative of the energy with respect to the variational parameter \( \alpha_i \) -We define also the derivative of the trial function (skipping the subindex \( T \)) as -

         
        -$$ -\bar{\Psi}_{i}=\frac{d\Psi}{d\alpha_i}. -$$ -

         
        -

        -
        - - -
        -

        Derivatives of the local energy

        -
        - -

        -The elements of the gradient of the local energy are then (using the chain rule and the hermiticity of the Hamiltonian) -

         
        -$$ -\bar{E}_{i}= 2\left( \langle \frac{\bar{\Psi}_{i}}{\Psi}E_L\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L \rangle\right). -$$ -

         
        - -From a computational point of view it means that you need to compute the expectation values of -

         
        -$$ -\langle \frac{\bar{\Psi}_{i}}{\Psi}E_L\rangle, -$$ -

         
        - -and -

         
        -$$ -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L\rangle -$$ -

         
        - -These integrals are evaluted using MC intergration (with all its possible error sources). -We can then use methods like stochastic gradient or other minimization methods to find the optimal variational parameters (I don't discuss this topic here, but these methods are very important in ML). -

        -
        - - -
        -

        How do we define our cost function?

        -
        - -

        -We have a model, our likelihood function. - -

        -How should we define the cost function? - - -

        -
        - - -
        -

        Meet the variance and its derivatives

        -
        -Why the variance? -

        -Suppose the trial function (our model) is the exact wave function. The action of the hamiltionan -on the wave function -

         
        -$$ - H\Psi = \mathrm{constant}\times \Psi, -$$ -

         
        - -The integral which defines various -expectation values involving moments of the hamiltonian becomes then -

         
        -$$ - \langle E^n \rangle = \langle H^n \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})H^n(\boldsymbol{R})\Psi(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})}= -\mathrm{constant}\times\frac{\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})}=\mathrm{constant}. -$$ -

         
        - -This gives an important information: If I want the variance, the exact wave function leads to zero variance! -The variance is defined as -

         
        -$$ -\sigma_E = \langle E^2\rangle - \langle E\rangle^2. -$$ -

         
        - -Variation is then performed by minimizing both the energy and the variance. - - -

        -
        - - -
        -

        The variance defines the cost function

        -
        - -

        -We can then take the derivatives of -

         
        -$$ -\sigma_E = \langle E^2\rangle - \langle E\rangle^2, -$$ -

         
        - -with respect to the variational parameters. The derivatives of the variance can then be used to defined the -so-called Hessian matrix, which in turn allows us to use minimization methods like Newton's method or -standard gradient methods. - -

        -This leads to however a more complicated expression, with obvious errors when evaluating integrals by Monte Carlo integration. Less used, see however Filippi and Umrigar. The expression becomes complicated -

         
        -$$ -\bar{E}_{ij} = 2\left[ \langle (\frac{\bar{\Psi}_{ij}}{\Psi}+\frac{\bar{\Psi}_{j}}{\Psi}\frac{\bar{\Psi}_{i}}{\Psi})(E_L-\langle E\rangle)\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\bar{E}_j-\langle \frac{\bar{\Psi}_{j}}{\Psi}\rangle\bar{E}_i\right] +\langle \frac{\bar{\Psi}_{i}}{\Psi}E_L{_j}\rangle +\langle \frac{\bar{\Psi}_{j}}{\Psi}E_L{_i}\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L{_j}\rangle \langle \frac{\bar{\Psi}_{j}}{\Psi}\rangle\langle E_L{_i}\rangle. -$$ -

         
        - -

        -Evaluating the cost function means having to evaluate the above second derivative of the energy. - - -

        -
        - - -
        -

        Boltzmann Machines

        +

        Boltzmann Machines

        Why use a generative model rather than the more well known discriminative deep neural networks (DNN)? @@ -1041,7 +218,7 @@ Why use a generative model rather than the more well known discriminative deep n

        -

        Some similarities and differences from DNNs

        +

        Some similarities and differences from DNNs

        1. Both use gradient-descent based learning procedures for minimizing cost functions
          2. @@ -1055,7 +232,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some
          -

          The structure of the RBM network

          +

          The structure of the RBM network





          @@ -1063,7 +240,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some
          -

          The network

          +

          The network

          The network layers: @@ -1076,7 +253,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some

          -

          Goals

          +

          Goals

          The goal of the hidden layer is to increase the model's expressive power. We 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). @@ -1102,13 +279,13 @@ Examples of this trick being employed in physics:

          -

          Joint distribution and the Energy function

          +

          Joint distribution and the Energy function

          The restricted Boltzmann machine is described by a Bolztmann distribution

           
          $$ \begin{align} P_{rbm}(\mathbf{x},\mathbf{h}) = \frac{1}{Z} e^{-\frac{1}{T_0}E(\mathbf{x},\mathbf{h})}, -\tag{2} +\tag{1} \end{align} $$

           
          @@ -1118,7 +295,7 @@ where \( Z \) is the normalization constant or partition function, defined as $$ \begin{align} Z = \int \int e^{-\frac{1}{T_0}E(\mathbf{x},\mathbf{h})} d\mathbf{x} d\mathbf{h}. -\tag{3} +\tag{2} \end{align} $$

           
          @@ -1128,7 +305,7 @@ It is common to ignore \( T_0 \) by setting it to one.

          -

          Network Elements

          +

          Network Elements

          The function \( E(\mathbf{x},\mathbf{h}) \) gives the energy of a @@ -1141,7 +318,7 @@ adjusting the energy function to best fit our problem.

          -

          Defining different types of RBMs

          +

          Defining different types of RBMs

          There are different variants of RBMs, and the differences lie in the types of visible and hidden units we choose as well as in the implementation of the energy function \( E(\mathbf{x},\mathbf{h}) \).

          @@ -1153,7 +330,7 @@ RBMs were first developed using binary units in both the visible and hidden laye $$ \begin{align} E(\mathbf{x}, \mathbf{h}) = - \sum_i^M x_i a_i- \sum_j^N b_j h_j - \sum_{i,j}^{M,N} x_i w_{ij} h_j, -\tag{4} +\tag{3} \end{align} $$

           
          @@ -1169,7 +346,7 @@ Another varient is the RBM where the visible units are Gaussian while the hidden $$ \begin{align} E(\mathbf{x}, \mathbf{h}) = \sum_i^M \frac{(x_i - a_i)^2}{2\sigma_i^2} - \sum_j^N b_j h_j - \sum_{i,j}^{M,N} \frac{x_i w_{ij} h_j}{\sigma_i^2}. -\tag{5} +\tag{4} \end{align} $$

           
          @@ -1178,7 +355,7 @@ $$

          -

          More about RBMs

          +

          More about RBMs

          1. Useful when we model continuous data (i.e., we wish \( \mathbf{x} \) to be continuous)
            2. @@ -1198,7 +375,7 @@ Other types of units include:
            -

            Sampling: Metropolis sampling

            +

            Sampling: Metropolis sampling

            In order to sample from the RBM probability distribution it is common to use Markov Chain Monte Carlo (MCMC) algorithms such as Metropolis-Hastings or Gibbs sampling.

            @@ -1207,7 +384,7 @@ Metropolis sampling starts by suggesting a new configuration \( \boldsymbol{x}^{ $$ \begin{align} A(\boldsymbol{x}^k \rightarrow \boldsymbol{x}^{k+1}) = \text{min} (1, \frac{P(\boldsymbol{x}^{k+1})}{P(\boldsymbol{x}^k)}), -\tag{6} +\tag{5} \end{align} $$

             
            @@ -1217,9 +394,9 @@ where we need the marginalized probability $$ \begin{align} P(\boldsymbol{x}) &= \sum_\mathbf{h} P_{rbm}(\mathbf{x}, \mathbf{h}) -\tag{7}\\ +\tag{6}\\ &= \frac{1}{Z}\sum_\mathbf{h} e^{-E(\mathbf{x}, \mathbf{h})}. -\tag{8} +\tag{7} \end{align} $$

             
            @@ -1227,7 +404,7 @@ $$

            -

            Sampling: Gibbs sampling

            +

            Sampling: Gibbs sampling

            In this method we sample from the joint probability \( P_{rbm} (\mathbf{x}, \mathbf{h}) \) by way of a two step sampling process. We alternately update the visible and hidden units. @@ -1239,9 +416,9 @@ The conditional probabilities are often referred to as the activitation function $$ \begin{align} P(h_j = 1 | \boldsymbol{x}) &= \frac{1}{1 + e^{-b_j - \sum_i x_i w_{ij}}} -\tag{9}\\ +\tag{8}\\ P(x_i = 1 | \boldsymbol{h}) &= \frac{1}{1 + e^{-a_j - \sum_j h_j w_{ij}}}, -\tag{10} +\tag{9} \end{align} $$

             
            @@ -1251,15 +428,15 @@ where we recognize the logistic sigmoid function \( \sigma (x) = 1/(1+exp(-x)) \

            -

            Gaussian RBM

            +

            Gaussian RBM

            For the Gaussian-Binary RBM the conditional probabilities are

             
            $$ \begin{align} P(x_i|\mathbf{h}) &= \mathcal{N}(x_i; a_i+ \sum_j h_j w_{ij}, \sigma^2) -\tag{11}\\ +\tag{10}\\ P(h_j=1|\mathbf{x}) &= \frac{1}{1+e^{-b_j-\frac{1}{\sigma^2} \sum_i x_i w_{ij}}}, -\tag{12} +\tag{11} \end{align} $$

             
            @@ -1269,7 +446,7 @@ while the visible units now follow a normal distribution, we see the hidden unit

            -

            Cost function

            +

            Cost function

            When working with a training dataset, the most common training approach is maximizing the log-likelihood of the training data. The log likelihood characterizes the log-probability of generating the observed data using our generative model. Using this method our cost function is chosen as the negative log-likelihood. The learning then consists of trying to find parameters that maximize the probability of the dataset, and is known as Maximum Likelihood Estimation (MLE). @@ -1278,9 +455,9 @@ Denoting the parameters as \( \boldsymbol{\theta} = a_1,...,a_M,b_1,...,b_N,w_{1 $$ \begin{align} \mathcal{L}(\{ \theta_i \}) &= \langle \text{log} P_\theta(\boldsymbol{x}) \rangle_{data} -\tag{13}\\ +\tag{12}\\ &= - \langle E(\boldsymbol{x}; \{ \theta_i\}) \rangle_{data} - \text{log} Z(\{ \theta_i\}), -\tag{14} +\tag{13} \end{align} $$

             
            @@ -1291,13 +468,13 @@ Our cost function is the negative log-likelihood, \( \mathcal{C}(\{ \theta_i \})

            -

            Optimization / Training

            +

            Optimization / Training

            The training procedure of choice often is Stochastic Gradient Descent (SGD). It consists of a series of iterations where we update the parameters according to the equation

             
            $$ \begin{align} \boldsymbol{\theta}_{k+1} = \boldsymbol{\theta}_k - \eta \nabla \mathcal{C} (\boldsymbol{\theta}_k) -\tag{15} +\tag{14} \end{align} $$

             
            @@ -1312,9 +489,9 @@ $$ \frac{\partial \mathcal{C}(\{ \theta_i\})}{\partial \theta_i} &= \langle \frac{\partial E(\boldsymbol{x}; \theta_i)}{\partial \theta_i} \rangle_{data} + \frac{\partial \text{log} Z(\{ \theta_i\})}{\partial \theta_i} -\tag{16}\\ +\tag{15}\\ &= \langle O_i(\boldsymbol{x}) \rangle_{data} - \langle O_i(\boldsymbol{x}) \rangle_{model}, -\tag{17} +\tag{16} \end{align} $$

             
            @@ -1324,7 +501,7 @@ where in order to simplify notation we defined the "operator" $$ \begin{align} O_i(\boldsymbol{x}) = \frac{\partial E(\boldsymbol{x}; \theta_i)}{\partial \theta_i}, -\tag{18} +\tag{17} \end{align} $$

             
            @@ -1334,7 +511,7 @@ and used the statistical mechanics relationship between expectation values and t $$ \begin{align} \langle O_i(\boldsymbol{x}) \rangle_{model} = \text{Tr} P_\theta(\boldsymbol{x})O_i(\boldsymbol{x}) = - \frac{\partial \text{log} Z(\{ \theta_i\})}{\partial \theta_i}. -\tag{19} +\tag{18} \end{align} $$

             
            @@ -1342,7 +519,7 @@ $$

            -

            More on RBMs

            +

            More on RBMs

            The data-dependent term in the gradient is known as the positive phase of the gradient, while the model-dependent term is known as the negative phase of the gradient. The aim of the training is to lower the energy of configurations that are near observed data points (increasing their probability), and raising the energy of configurations that are far from observed data points (decreasing their probability). @@ -1353,12 +530,12 @@ The gradient of the negative log-likelihood cost function of a Binary-Binary RBM $$ \begin{align} \frac{\partial \mathcal{C} (w_{ij}, a_i, b_j)}{\partial w_{ij}} =& \langle x_i h_j \rangle_{data} - \langle x_i h_j \rangle_{model} -\tag{20}\\ +\tag{19}\\ \frac{\partial \mathcal{C} (w_{ij}, a_i, b_j)}{\partial a_{ij}} =& \langle x_i \rangle_{data} - \langle x_i \rangle_{model} -\tag{21}\\ +\tag{20}\\ \frac{\partial \mathcal{C} (w_{ij}, a_i, b_j)}{\partial b_{ij}} =& \langle h_i \rangle_{data} - \langle h_i \rangle_{model}. -\tag{22}\\ -\tag{23} +\tag{21}\\ +\tag{22} \end{align} $$

             
            @@ -1368,7 +545,7 @@ To get the expecation values with respect to the data, we set the visib

            -

            Which sampling to use

            +

            Which sampling to use

            To get the expectation values with respect to the model, we use Gibbs sampling. We can either initialize the \( \boldsymbol{x} \) randomly or with a training sample. While we ideally want a large number of Gibbs iterations \( n\rightarrow n \), one might decide to truncate it earlier for efficiency. Doing this while having intialized \( \boldsymbol{x} \) with a training data vector is referred to as contrastive divergence (CD), because one is then closer to approximating the gradient of this function than the negative log-likelihood. The contrastive divergence function is the difference between two Kullback-Leibler divergences (also called relative entropy), which measure how one probability distribution diverges from a second, expected probability distribution (in this case the estimated one from the ground truth one). @@ -1376,7 +553,7 @@ To get the expectation values with respect to the model, we use Gibbs s

            -

            RBMs for the quantum many body problem

            +

            RBMs for the quantum many body problem

            The idea of applying RBMs to quantum many body problems was presented by G. Carleo and M. Troyer, working with ETH Zurich and Microsoft Research. @@ -1400,7 +577,7 @@ Some of their motivation included

            -

            Choose the right RBM

            +

            Choose the right RBM

            Carleo and Troyer applied the RBM to the quantum mechanical spin lattice systems of the Ising model and Heisenberg model, with encouraging results. Our goal is to test the method on systems of moving particles. For the spin lattice systems it was natural to use a binary-binary RBM, with the nodes taking values of 1 and -1. For moving particles, on the other hand, we want the visible nodes to be continuous, representing position coordinates. Thus, we start by choosing a Gaussian-binary RBM, where the visible nodes are continuous and hidden nodes take on values of 0 or 1. If eventually we would like the hidden nodes to be continuous as well the rectified linear units seem like the most relevant choice. @@ -1408,13 +585,13 @@ Carleo and Troyer applied the RBM to the quantum mechanical spin lattice systems

            -

            Representing the wave function

            +

            Representing the wave function

            The wavefunction should be a probability amplitude depending on \( \boldsymbol{x} \). The RBM model is given by the joint distribution of \( \boldsymbol{x} \) and \( \boldsymbol{h} \)

             
            $$ \begin{align} F_{rbm}(\mathbf{x},\mathbf{h}) = \frac{1}{Z} e^{-\frac{1}{T_0}E(\mathbf{x},\mathbf{h})}. -\tag{24} +\tag{23} \end{align} $$

             
            @@ -1424,9 +601,9 @@ To find the marginal distribution of \( \boldsymbol{x} \) we set: $$ \begin{align} F_{rbm}(\mathbf{x}) &= \sum_\mathbf{h} F_{rbm}(\mathbf{x}, \mathbf{h}) -\tag{25}\\ +\tag{24}\\ &= \frac{1}{Z}\sum_\mathbf{h} e^{-E(\mathbf{x}, \mathbf{h})}. -\tag{26} +\tag{25} \end{align} $$

             
            @@ -1436,14 +613,14 @@ Now this is what we use to represent the wave function, calling it a neural-netw $$ \begin{align} \Psi (\mathbf{X}) &= F_{rbm}(\mathbf{x}) -\tag{27}\\ +\tag{26}\\ &= \frac{1}{Z}\sum_{\boldsymbol{h}} e^{-E(\mathbf{x}, \mathbf{h})} -\tag{28}\\ +\tag{27}\\ &= \frac{1}{Z} \sum_{\{h_j\}} e^{-\sum_i^M \frac{(x_i - a_i)^2}{2\sigma^2} + \sum_j^N b_j h_j + \sum_{i,j}^{M,N} \frac{x_i w_{ij} h_j}{\sigma^2}} -\tag{29}\\ +\tag{28}\\ &= \frac{1}{Z} e^{-\sum_i^M \frac{(x_i - a_i)^2}{2\sigma^2}} \prod_j^N (1 + e^{b_j + \sum_i^M \frac{x_i w_{ij}}{\sigma^2}}). -\tag{30}\\ -\tag{31} +\tag{29}\\ +\tag{30} \end{align} $$

             
            @@ -1451,14 +628,14 @@ $$

            -

            Choose the cost function

            +

            Choose the cost function

            Now we don't necessarily have training data (unless we generate it by using some other method). However, what we do have is the variational principle which allows us to obtain the ground state wave function by minimizing the expectation value of the energy of a trial wavefunction (corresponding to the untrained NQS). Similarly to the traditional variational Monte Carlo method then, it is the local energy we wish to minimize. The gradient to use for the stochastic gradient descent procedure is

             
            $$ \begin{align} G_i = \frac{\partial \langle E_L \rangle}{\partial \theta_i} = 2(\langle E_L \frac{1}{\Psi}\frac{\partial \Psi}{\partial \theta_i} \rangle - \langle E_L \rangle \langle \frac{1}{\Psi}\frac{\partial \Psi}{\partial \theta_i} \rangle ), -\tag{32} +\tag{31} \end{align} $$

             
            @@ -1468,7 +645,7 @@ where the local energy is given by $$ \begin{align} E_L = \frac{1}{\Psi} \hat{\mathbf{H}} \Psi. -\tag{33} +\tag{32} \end{align} $$

             
            @@ -1476,7 +653,7 @@ $$

            -

            Running the codes

            +

            Running the codes

            @@ -1494,7 +671,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat

            -

            Energy as function of iterations, \( N=2 \) electrons

            +

            Energy as function of iterations, \( N=2 \) electrons

            @@ -1504,7 +681,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat

            -

            Energy as function of iterations, \( N=6 \) electrons

            +

            Energy as function of iterations, \( N=6 \) electrons

            @@ -1514,7 +691,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat

            -

            Conclusions and where do we stand

            +

            Conclusions and where do we stand

              diff --git a/doc/pub/BM/html/BM-solarized.html b/doc/pub/BM/html/BM-solarized.html index cacbfbcfb..f2642008e 100644 --- a/doc/pub/BM/html/BM-solarized.html +++ b/doc/pub/BM/html/BM-solarized.html @@ -6,9 +6,9 @@ Automatically generated HTML file from DocOnce source - + -Solving quantum mechanical problems with Machine Learning and Boltzmann machines +Machine Learning and Boltzmann machines @@ -60,95 +60,44 @@ div { text-align: justify; text-justify: inter-word; } @@ -174,7 +123,7 @@ MathJax.Hub.Config({ -

              Solving quantum mechanical problems with Machine Learning and Boltzmann machines

              +

              Machine Learning and Boltzmann machines

              @@ -189,23 +138,15 @@ MathJax.Hub.Config({

              Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University and Department of Physics, University of Oslo, Norway

              -

              Jul 26, 2018

              +

              Nov 19, 2018












              -

              What is this talk about?

              +

              Unsupervised learning, ovrarching aims

              -The main aim is to give you a short and pedestrian introduction to how we can use Machine Learning methods -to solve quantum mechanical many-body problems. And why this could be of interest. - -

              -The hope is that after this talk you have gotten the basic ideas to get you started. Peeping into https://github.com/mhjensenseminars/MachineLearningTalk, you'll find a Jupyter notebook, slides, codes etc that will allow you to reproduce the simulations discussed here, and perhaps run your own very first calculations. - -

              -Furthermore, since we will also talk about Monte Carlo simulations, the exciting idea of Monte Carlo optimization and Machine Learning can be used to mesh strategies for reducing the curse of dimensionality.

              @@ -214,188 +155,7 @@ Furthermore, since we will also talk about Monte Carlo simulations, the exciting











              -

              Why?

              - -

              -How can we avoid the dimensionality curse? Many possibilities - -

                -
              1. smarter basis functions
                2. -
              2. resummation of specific correlations
                3. -
              3. stochastic sampling of high-lying states (stochastic FCI, CC and soon SRG)
                4. -
              4. many more
                5. -
              - -Machine Learning and Quantum Computing hold also great promise in tackling the -ever increasing dimensionalities. Here we will focus on Machine Learning. - -

              -









              - -

              Overview

              -
              - -

              - -

                -
              • Short intro to Machine Learning
                • -
              • Variational Monte Carlo (Markov Chain Monte Carlo, \( \mathrm{MC}^2 \)) and two-electron quantum dots, solving quantum mechanical problems in a stochastic way. It will serve as our motivation for switching to Machine Learning.
                • -
              • From Variational Monte Carlo to Boltzmann Machines and Machine Learning
                • -
              -
              - - -

              -









              - -

              What are the Machine Learning calculations here based on?

              - -

              -This work is inspired by the idea of representing the wave function with -a restricted Boltzmann machine (RBM), presented recently by G. Carleo and M. Troyer, Science 355, Issue 6325, pp. 602-606 (2017). They -named such a wave function/network a neural network quantum state (NQS). In their article they apply it to the quantum mechanical -spin lattice systems of the Ising model and Heisenberg model, with -encouraging results. - -

              -Thanks to Vilde Flugsrud and Alfred Alocias Mariadason for many discussions and interpretations of the above paper. - -

              -









              - -

              A new world

              -
              - -

              -Machine learning (ML) 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 ML are - -

              - -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. -
              - - -

              -









              - -

              Lots of room for creativity

              -
              - -

              -Not all the -algorithms and methods can be given a rigorous mathematical -justification, opening up thereby for experimenting -and trial and error and thereby exciting new developments. -

              - - -

              -

              - -

              -A solid command of linear algebra, multivariate theory, -probability theory, statistical data analysis, -understanding errors and Monte Carlo methods is important in order to understand many of the -various algorithms and methods. -

              - - -

              -

              -Job market, a personal statement. -

              -A familiarity with ML is almost becoming a prerequisite for many of the most exciting employment -opportunities. Add quantum computing and there you are. -

              - - -

              - - -

              Some members of the ML family

              - -

              -

              - -

              - -

                -
              1. Linear regression and its variants, in essence polynomial regression
                2. -
              2. Decision tree algorithms, from simpler to more complex ones like random forests
                3. -
              3. Nearest neighbors models
                4. -
              4. Bayesian statistics and regression
                5. -
              5. Support vector machines and finally various variants of
                6. -
              6. Artifical neural networks and deep learning
                7. -
              7. ....
                8. -
              -
              - - -

              -









              - -

              What are the basic ingredients?

              -
              - -

              -Almost every problem in ML and data science starts with the same ingredients: - -

                -
              • The dataset \( \mathbf{x} \) (could be some observable quantity of the system we are studying)
                • -
              • A model which is a function of a set of parameters \( \mathbf{\alpha} \) that relates to the dataset, say a likelihood function \( p(\mathbf{x}\vert \mathbf{\alpha}) \) or just a simple model \( f(\mathbf{\alpha}) \)
                • -
              • A so-called cost function \( \mathcal{C} (\mathbf{x}, f(\mathbf{\alpha})) \) which allows us to decide how well our model represents the dataset.
                • -
              - -We seek to minimize the function \( \mathcal{C} (\mathbf{x}, f(\mathbf{\alpha})) \) by finding the parameter values which minimize \( \mathcal{C} \). This leads to various minimization algorithms. -
              - - -

              -









              - -

              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 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. - -

              -You will notice however that many of the basic ideas discussed do come from Physics! -

              - - -

              -









              - -

              Types of Machine Learning

              +

              Types of Machine Learning

              @@ -424,615 +184,7 @@ Some of the most common tasks are:











              -

              References

              -
              - -

              - -

              -
              - - -

              -









              - -

              Another interesting article

              - -

              -Here we will use so-called reduced Boltzmann Machines to simulate quantum many-body problems. For Monte Carlo aficionados, there is a very close similarity with what are called shadow wave functions, see the work of Pederiva and Kalos and collaborators, Phys Rev. E 90, 053304 (2014). - -

              -









              - -

              Just for fun illustration of scikit-learn: Hudson bay data

              -
              - -

              -Most mammalian predators rely on a variety of prey, which complicates mathematical modeling; however, a few predators have become highly specialized and seek almost exclusively a single prey species. An example of this simplified predator-prey interaction is seen in Canadian northern forests, where the populations of the lynx and the snowshoe hare are intertwined in a life and death struggle. - -

              -One reason that this particular system has been so extensively studied is that the Hudson Bay company kept careful records of all furs from the early 1800s into the 1900s. The records for the furs collected by the Hudson Bay company showed distinct oscillations (approximately 12 year periods), suggesting that these species caused almost periodic fluctuations of each other's populations. The table here shows data from 1900 to 1920. - - -

              - - -

              -









              - -

              Hares and lynx in Hudson bay from 1900 to 1920

              - -

              -



              - -

              -









              - -

              Linear regression in Python

              -
              - -

              -

              - - -

              import numpy as np
-import matplotlib.pyplot as plt
-from IPython.display import display
-import sklearn
-from sklearn.linear_model import LinearRegression
-from sklearn.tree import DecisionTreeRegressor
-
-
-data = np.loadtxt('src/Hudson_Bay.csv', delimiter=',', skiprows=1)
-x = data[:,0]
-y = data[:,1]
-line = np.linspace(1900,1920,1000,endpoint=False).reshape(-1,1)
-reg = DecisionTreeRegressor(min_samples_split=3).fit(x.reshape(-1,1),y.reshape(-1,1))
-plt.plot(line, reg.predict(line), label="decision tree")
-regline = LinearRegression().fit(x.reshape(-1,1),y.reshape(-1,1))
-plt.plot(line, regline.predict(line), label= "Linear Regression")
-plt.plot(x, y, label= "Linear Regression")
-plt.show()
-
              - -
              - - -

              -









              - -

              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. - -

              -This behaviour has inspired a simple mathematical model for an artificial neuron. - -$$ - y = f\left(\sum_{i=1}^n w_ix_i\right) = f(u) -$$ - -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. -

              - - -

              -









              - -

              A simple perceptron model

              - -

              -



              - -

              -









              - -

              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. - -

              -



              -
              - - -

              -









              - -

              The system: two electrons in a harmonic oscillator trap in two dimensions

              - -

              -The Hamiltonian of the quantum dot is given by -$$ \hat{H} = \hat{H}_0 + \hat{V}, -$$ - -where \( \hat{H}_0 \) is the many-body HO Hamiltonian, and \( \hat{V} \) is the -inter-electron Coulomb interactions. In dimensionless units, -$$ \hat{V}= \sum_{i < j}^N \frac{1}{r_{ij}}, -$$ - -with \( r_{ij}=\sqrt{\mathbf{r}_i^2 - \mathbf{r}_j^2} \). - -

              -This leads to the separable Hamiltonian, with the relative motion part given by (\( r_{ij}=r \)) -$$ -\hat{H}_r=-\nabla^2_r + \frac{1}{4}\omega^2r^2+ \frac{1}{r}, -$$ - -plus a standard Harmonic Oscillator problem for the center-of-mass motion. -This system has analytical solutions in two and three dimensions (M. Taut 1993 and 1994). - -

              -









              - -

              Quantum Monte Carlo Motivation

              -
              - -

              -Given a hamiltonian \( H \) and a trial wave function \( \Psi_T \), the variational principle states that the expectation value of \( \langle H \rangle \), defined through -$$ - \langle E \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})H(\boldsymbol{R})\Psi_T(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}, -$$ - -is an upper bound to the ground state energy \( E_0 \) of the hamiltonian \( H \), that is -$$ - E_0 \le \langle H \rangle . -$$ - -In general, the integrals involved in the calculation of various expectation values are multi-dimensional ones. Traditional integration methods such as the Gauss-Legendre will not be adequate for say the computation of the energy of a many-body system. -

              - - -

              -









              - -

              Quantum Monte Carlo Motivation

              -
              -Basic steps. -

              -Choose a trial wave function -\( \psi_T(\boldsymbol{R}) \). -$$ - P(\boldsymbol{R},\boldsymbol{\alpha})= \frac{\left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2}{\int \left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2d\boldsymbol{R}}. -$$ - -This is our model, or likelihood/probability distribution function (PDF). It depends on some variational parameters \( \boldsymbol{\alpha} \). -The approximation to the expectation value of the Hamiltonian is now -$$ - \langle E[\boldsymbol{\alpha}] \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R},\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_T(\boldsymbol{R},\boldsymbol{\alpha})} - {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R},\boldsymbol{\alpha})\Psi_T(\boldsymbol{R},\boldsymbol{\alpha})}. -$$ -

              - - -

              -









              - -

              Quantum Monte Carlo Motivation

              -
              -Define a new quantity. -

              -$$ - E_L(\boldsymbol{R},\boldsymbol{\alpha})=\frac{1}{\psi_T(\boldsymbol{R},\boldsymbol{\alpha})}H\psi_T(\boldsymbol{R},\boldsymbol{\alpha}), -$$ - -called the local energy, which, together with our trial PDF yields -$$ - E[\boldsymbol{\alpha}]=\int P(\boldsymbol{R})E_L(\boldsymbol{R},\boldsymbol{\alpha}) d\boldsymbol{R}\approx \frac{1}{N}\sum_{i=1}^NE_L(\boldsymbol{R_i},\boldsymbol{\alpha}) -$$ - -with \( N \) being the number of Monte Carlo samples. -

              - - -

              -









              - -

              Quantum Monte Carlo

              -
              - -

              -The Algorithm for performing a variational Monte Carlo calculations runs thus as this - -

                -
              • Initialisation: Fix the number of Monte Carlo steps. Choose an initial \( \boldsymbol{R} \) and variational parameters \( \alpha \) and calculate \( \left|\psi_T(\boldsymbol{R},\boldsymbol{\alpha})\right|^2 \).
                • -
              • Initialise the energy and the variance and start the Monte Carlo calculation by looping over trials.
                • - -
                  -
                • Calculate a trial position \( \boldsymbol{R}_p=\boldsymbol{R}+r*step \) where \( r \) is a random variable \( r \in [0,1] \).
                  • -
                • Metropolis algorithm to accept or reject this move \( w = P(\boldsymbol{R}_p,\boldsymbol{\alpha})/P(\boldsymbol{R},\boldsymbol{\alpha}) \).
                  • -
                • If the step is accepted, then we set \( \boldsymbol{R}=\boldsymbol{R}_p \).
                  • -
                • Update averages
                  • -
                - -
              • Finish and compute final averages.
                • -
              - -Observe that the jumping in space is governed by the variable step. This is often called brute-force sampling. -Need importance sampling to get more relevant sampling. -
              - - -

              -









              - -

              The trial wave function

              -
              - -

              -We want to perform a Variational Monte Carlo calculation of the ground state of two electrons in a quantum dot well with different oscillator energies, assuming total spin \( S=0 \). -Our trial wave function has the following form -$$ -\begin{equation} - \psi_{T}(\boldsymbol{r}_1,\boldsymbol{r}_2) = - C\exp{\left(-\alpha_1\omega(r_1^2+r_2^2)/2\right)} - \exp{\left(\frac{r_{12}}{(1+\alpha_2 r_{12})}\right)}, -\label{eq:trial} -\end{equation} -$$ - -where the $\alpha$s represent our variational parameters, two in this case. - -

              -Why does the trial function look like this? How did we get there? This will be our main motivation for switching to -Machine Learning. - - -

              - - -

              -









              - -

              The correlation part of the wave function

              - -

              -To find an ansatz for the correlated part of the wave function, it is useful to rewrite the two-particle -local energy in terms of the relative and center-of-mass motion. -Let us denote the distance between the two electrons as -\( r_{12} \). We omit the center-of-mass motion since we are only interested in the case when -\( r_{12} \rightarrow 0 \). The contribution from the center-of-mass (CoM) variable \( {\bf R}_{\mathrm{CoM}} \) -gives only a finite contribution. -We focus only on the terms that are relevant for \( r_{12} \) and for three dimensions. The relevant local energy becomes then -$$ -\lim_{r_{12} \rightarrow 0}E_L(R)= - \frac{1}{{\cal R}_T(r_{12})}\left(2\frac{d^2}{dr_{ij}^2}+\frac{4}{r_{ij}}\frac{d}{dr_{ij}}+ -\frac{2}{r_{ij}}-\frac{l(l+1)}{r_{ij}^2}+2E -\right){\cal R}_T(r_{12}) = 0. -$$ - -Set \( l=0 \) and we have the so-called cusp condition -$$ -\frac{d {\cal R}_T(r_{12})}{dr_{12}} = -\frac{1}{2(l+1)} -{\cal R}_T(r_{12})\qquad r_{12}\to 0 -$$ - -

              -









              - -

              Resulting ansatz

              -The above results in -$$ -{\cal R}_T \propto \exp{(r_{ij}/2)}, -$$ - -for anti-parallel spins and -$$ -{\cal R}_T \propto \exp{(r_{ij}/4)}, -$$ - -for anti-parallel spins. -This is the so-called cusp condition for the relative motion, resulting in a minimal requirement -for the correlation part of the wave fuction. -For general systems containing more than say two electrons, we have this -condition for each electron pair \( ij \). - -

              -









              - -

              The VMC code

              - -

              - - -

              # Importing various packages
-from math import exp, sqrt
-from random import random, seed
-import numpy as np
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-from matplotlib import cm
-from matplotlib.ticker import LinearLocator, FormatStrFormatter
-import sys
-
-#Trial wave function for quantum dots in two dims
-def WaveFunction(r,alpha,beta):
-    r1 = r[0,0]**2 + r[0,1]**2
-    r2 = r[1,0]**2 + r[1,1]**2
-    r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
-    deno = r12/(1+beta*r12)
-    return exp(-0.5*alpha*(r1+r2)+deno)
-
-#Local energy  for quantum dots in two dims, using analytical local energy
-def LocalEnergy(r,alpha,beta):
-    
-    r1 = (r[0,0]**2 + r[0,1]**2)
-    r2 = (r[1,0]**2 + r[1,1]**2)
-    r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
-    deno = 1.0/(1+beta*r12)
-    deno2 = deno*deno
-    return 0.5*(1-alpha*alpha)*(r1 + r2) +2.0*alpha + 1.0/r12+deno2*(alpha*r12-deno2+2*beta*deno-1.0/r12)
-
-# The Monte Carlo sampling with the Metropolis algo
-def MonteCarloSampling():
-
-    NumberMCcycles= 100000
-    StepSize = 1.0
-    # positions
-    PositionOld = np.zeros((NumberParticles,Dimension), np.double)
-    PositionNew = np.zeros((NumberParticles,Dimension), np.double)
-    # seed for rng generator
-    seed()
-    # start variational parameter
-    alpha = 0.9
-    for ia in range(MaxVariations):
-        alpha += .025
-        AlphaValues[ia] = alpha
-        beta = 0.2 
-        for jb in range(MaxVariations):
-            beta += .01
-            BetaValues[jb] = beta
-            energy = energy2 = 0.0
-            DeltaE = 0.0
-            #Initial position
-            for i in range(NumberParticles):
-                for j in range(Dimension):
-                    PositionOld[i,j] = StepSize * (random() - .5)
-            wfold = WaveFunction(PositionOld,alpha,beta)
-
-            #Loop over MC MCcycles
-            for MCcycle in range(NumberMCcycles):
-                #Trial position
-                for i in range(NumberParticles):
-                    for j in range(Dimension):
-                        PositionNew[i,j] = PositionOld[i,j] + StepSize * (random() - .5)
-                wfnew = WaveFunction(PositionNew,alpha,beta)
-
-                #Metropolis test to see whether we accept the move
-                if random() < wfnew**2 / wfold**2:
-                   PositionOld = PositionNew.copy()
-                   wfold = wfnew
-                   DeltaE = LocalEnergy(PositionOld,alpha,beta)
-                energy += DeltaE
-                energy2 += DeltaE**2
-
-            #We calculate mean, variance and error ...
-            energy /= NumberMCcycles
-            energy2 /= NumberMCcycles
-            variance = energy2 - energy**2
-            error = sqrt(variance/NumberMCcycles)
-            Energies[ia,jb] = energy    
-    return Energies, AlphaValues, BetaValues
-
-
-#Here starts the main program with variable declarations
-NumberParticles = 2
-Dimension = 2
-MaxVariations = 10
-Energies = np.zeros((MaxVariations,MaxVariations))
-AlphaValues = np.zeros(MaxVariations)
-BetaValues = np.zeros(MaxVariations)
-(Energies, AlphaValues, BetaValues) = MonteCarloSampling()
-
-# Prepare for plots
-fig = plt.figure()
-ax = fig.gca(projection='3d')
-# Plot the surface.
-X, Y = np.meshgrid(AlphaValues, BetaValues)
-surf = ax.plot_surface(X, Y, Energies,cmap=cm.coolwarm,linewidth=0, antialiased=False)
-# Customize the z axis.
-zmin = np.matrix(Energies).min()
-zmax = np.matrix(Energies).max()
-ax.set_zlim(zmin, zmax)
-ax.set_xlabel(r'$\alpha$')
-ax.set_ylabel(r'$\beta$')
-ax.set_zlabel(r'$\langle E \rangle$')
-ax.zaxis.set_major_locator(LinearLocator(10))
-ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
-# Add a color bar which maps values to colors.
-fig.colorbar(surf, shrink=0.5, aspect=5)
-plt.show()
-
              -

              -









              - -

              Technical aspect, improvements and how to define the cost function

              -
              - -

              - -

              -The above procedure is not the smartest one. Looping over all variational parameters becomes expensive. -Also, we don't use importance sampling and optimizations of the standard deviation (blocking, bootstrap, jackknife). -Such codes are included in the above Github address. - -

              -We can also be smarter and use minimization methods to find the optimal variational parameters with fewer Monte Carlo cycles and then -fire up our heavy artillery. - -

              -One way to achieve this is to minimize the energy as function of the variational parameters. - - -

              - - -

              -









              - -

              Energy derivatives

              -
              - -

              -To find the derivatives of the local energy expectation value as function of the variational parameters, we can use the chain rule and the hermiticity of the Hamiltonian. - -

              -Let us define -$$ -\bar{E}_{\alpha_i}=\frac{d\langle E_L\rangle}{d\alpha_i}. -$$ - -as the derivative of the energy with respect to the variational parameter \( \alpha_i \) -We define also the derivative of the trial function (skipping the subindex \( T \)) as -$$ -\bar{\Psi}_{i}=\frac{d\Psi}{d\alpha_i}. -$$ -

              - - -

              -









              - -

              Derivatives of the local energy

              -
              - -

              -The elements of the gradient of the local energy are then (using the chain rule and the hermiticity of the Hamiltonian) -$$ -\bar{E}_{i}= 2\left( \langle \frac{\bar{\Psi}_{i}}{\Psi}E_L\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L \rangle\right). -$$ - -From a computational point of view it means that you need to compute the expectation values of -$$ -\langle \frac{\bar{\Psi}_{i}}{\Psi}E_L\rangle, -$$ - -and -$$ -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L\rangle -$$ - -These integrals are evaluted using MC intergration (with all its possible error sources). -We can then use methods like stochastic gradient or other minimization methods to find the optimal variational parameters (I don't discuss this topic here, but these methods are very important in ML). -

              - - -

              -









              - -

              How do we define our cost function?

              -
              - -

              -We have a model, our likelihood function. - -

              -How should we define the cost function? - - -

              - - -

              -









              - -

              Meet the variance and its derivatives

              -
              -Why the variance? -

              -Suppose the trial function (our model) is the exact wave function. The action of the hamiltionan -on the wave function -$$ - H\Psi = \mathrm{constant}\times \Psi, -$$ - -The integral which defines various -expectation values involving moments of the hamiltonian becomes then -$$ - \langle E^n \rangle = \langle H^n \rangle = - \frac{\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})H^n(\boldsymbol{R})\Psi(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})}= -\mathrm{constant}\times\frac{\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})} - {\int d\boldsymbol{R}\Psi^{\ast}(\boldsymbol{R})\Psi(\boldsymbol{R})}=\mathrm{constant}. -$$ - -This gives an important information: If I want the variance, the exact wave function leads to zero variance! -The variance is defined as -$$ -\sigma_E = \langle E^2\rangle - \langle E\rangle^2. -$$ - -Variation is then performed by minimizing both the energy and the variance. - - -

              - - -

              -









              - -

              The variance defines the cost function

              -
              - -

              -We can then take the derivatives of -$$ -\sigma_E = \langle E^2\rangle - \langle E\rangle^2, -$$ - -with respect to the variational parameters. The derivatives of the variance can then be used to defined the -so-called Hessian matrix, which in turn allows us to use minimization methods like Newton's method or -standard gradient methods. - -

              -This leads to however a more complicated expression, with obvious errors when evaluating integrals by Monte Carlo integration. Less used, see however Filippi and Umrigar. The expression becomes complicated -$$ -\bar{E}_{ij} = 2\left[ \langle (\frac{\bar{\Psi}_{ij}}{\Psi}+\frac{\bar{\Psi}_{j}}{\Psi}\frac{\bar{\Psi}_{i}}{\Psi})(E_L-\langle E\rangle)\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\bar{E}_j-\langle \frac{\bar{\Psi}_{j}}{\Psi}\rangle\bar{E}_i\right] +\langle \frac{\bar{\Psi}_{i}}{\Psi}E_L{_j}\rangle +\langle \frac{\bar{\Psi}_{j}}{\Psi}E_L{_i}\rangle -\langle \frac{\bar{\Psi}_{i}}{\Psi}\rangle\langle E_L{_j}\rangle \langle \frac{\bar{\Psi}_{j}}{\Psi}\rangle\langle E_L{_i}\rangle. -$$ - -

              -Evaluating the cost function means having to evaluate the above second derivative of the energy. - - -

              - - -

              -









              - -

              Boltzmann Machines

              +

              Boltzmann Machines

              Why use a generative model rather than the more well known discriminative deep neural networks (DNN)? @@ -1051,7 +203,7 @@ Why use a generative model rather than the more well known discriminative deep n









              -

              Some similarities and differences from DNNs

              +

              Some similarities and differences from DNNs

              1. Both use gradient-descent based learning procedures for minimizing cost functions
                2. @@ -1064,7 +216,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some











                -

                The structure of the RBM network

                +

                The structure of the RBM network





                @@ -1072,7 +224,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some











                -

                The network

                +

                The network

                The network layers: @@ -1084,7 +236,7 @@ History: The RBM was developed by amongst others Geoffrey Hinton, called by some









                -

                Goals

                +

                Goals

                The goal of the hidden layer is to increase the model's expressive power. We 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). @@ -1108,7 +260,7 @@ Examples of this trick being employed in physics:









                -

                Joint distribution and the Energy function

                +

                Joint distribution and the Energy function

                The restricted Boltzmann machine is described by a Bolztmann distribution $$ \begin{align} @@ -1130,7 +282,7 @@ It is common to ignore \( T_0 \) by setting it to one.











                -

                Network Elements

                +

                Network Elements

                The function \( E(\mathbf{x},\mathbf{h}) \) gives the energy of a @@ -1143,7 +295,7 @@ adjusting the energy function to best fit our problem.











                -

                Defining different types of RBMs

                +

                Defining different types of RBMs

                There are different variants of RBMs, and the differences lie in the types of visible and hidden units we choose as well as in the implementation of the energy function \( E(\mathbf{x},\mathbf{h}) \).

                @@ -1181,7 +333,7 @@ $$











                -

                More about RBMs

                +

                More about RBMs

                1. Useful when we model continuous data (i.e., we wish \( \mathbf{x} \) to be continuous)
                  2. @@ -1199,7 +351,7 @@ Other types of units include:









                  -

                  Sampling: Metropolis sampling

                  +

                  Sampling: Metropolis sampling

                  In order to sample from the RBM probability distribution it is common to use Markov Chain Monte Carlo (MCMC) algorithms such as Metropolis-Hastings or Gibbs sampling.

                  @@ -1224,7 +376,7 @@ $$











                  -

                  Sampling: Gibbs sampling

                  +

                  Sampling: Gibbs sampling

                  In this method we sample from the joint probability \( P_{rbm} (\mathbf{x}, \mathbf{h}) \) by way of a two step sampling process. We alternately update the visible and hidden units. @@ -1246,7 +398,7 @@ where we recognize the logistic sigmoid function \( \sigma (x) = 1/(1+exp(-x)) \











                  -

                  Gaussian RBM

                  +

                  Gaussian RBM

                  For the Gaussian-Binary RBM the conditional probabilities are $$ \begin{align} @@ -1262,7 +414,7 @@ while the visible units now follow a normal distribution, we see the hidden unit











                  -

                  Cost function

                  +

                  Cost function

                  When working with a training dataset, the most common training approach is maximizing the log-likelihood of the training data. The log likelihood characterizes the log-probability of generating the observed data using our generative model. Using this method our cost function is chosen as the negative log-likelihood. The learning then consists of trying to find parameters that maximize the probability of the dataset, and is known as Maximum Likelihood Estimation (MLE). @@ -1282,7 +434,7 @@ Our cost function is the negative log-likelihood, \( \mathcal{C}(\{ \theta_i \})











                  -

                  Optimization / Training

                  +

                  Optimization / Training

                  The training procedure of choice often is Stochastic Gradient Descent (SGD). It consists of a series of iterations where we update the parameters according to the equation $$ \begin{align} @@ -1325,7 +477,7 @@ $$

                  -

                  More on RBMs

                  +

                  More on RBMs

                  The data-dependent term in the gradient is known as the positive phase of the gradient, while the model-dependent term is known as the negative phase of the gradient. The aim of the training is to lower the energy of configurations that are near observed data points (increasing their probability), and raising the energy of configurations that are far from observed data points (decreasing their probability). @@ -1349,7 +501,7 @@ To get the expecation values with respect to the data, we set the visib











                  -

                  Which sampling to use

                  +

                  Which sampling to use

                  To get the expectation values with respect to the model, we use Gibbs sampling. We can either initialize the \( \boldsymbol{x} \) randomly or with a training sample. While we ideally want a large number of Gibbs iterations \( n\rightarrow n \), one might decide to truncate it earlier for efficiency. Doing this while having intialized \( \boldsymbol{x} \) with a training data vector is referred to as contrastive divergence (CD), because one is then closer to approximating the gradient of this function than the negative log-likelihood. The contrastive divergence function is the difference between two Kullback-Leibler divergences (also called relative entropy), which measure how one probability distribution diverges from a second, expected probability distribution (in this case the estimated one from the ground truth one). @@ -1357,7 +509,7 @@ To get the expectation values with respect to the model, we use Gibbs s











                  -

                  RBMs for the quantum many body problem

                  +

                  RBMs for the quantum many body problem

                  The idea of applying RBMs to quantum many body problems was presented by G. Carleo and M. Troyer, working with ETH Zurich and Microsoft Research. @@ -1381,7 +533,7 @@ Some of their motivation included









                  -

                  Choose the right RBM

                  +

                  Choose the right RBM

                  Carleo and Troyer applied the RBM to the quantum mechanical spin lattice systems of the Ising model and Heisenberg model, with encouraging results. Our goal is to test the method on systems of moving particles. For the spin lattice systems it was natural to use a binary-binary RBM, with the nodes taking values of 1 and -1. For moving particles, on the other hand, we want the visible nodes to be continuous, representing position coordinates. Thus, we start by choosing a Gaussian-binary RBM, where the visible nodes are continuous and hidden nodes take on values of 0 or 1. If eventually we would like the hidden nodes to be continuous as well the rectified linear units seem like the most relevant choice. @@ -1389,7 +541,7 @@ Carleo and Troyer applied the RBM to the quantum mechanical spin lattice systems











                  -

                  Representing the wave function

                  +

                  Representing the wave function

                  The wavefunction should be a probability amplitude depending on \( \boldsymbol{x} \). The RBM model is given by the joint distribution of \( \boldsymbol{x} \) and \( \boldsymbol{h} \) $$ \begin{align} @@ -1426,7 +578,7 @@ $$











                  -

                  Choose the cost function

                  +

                  Choose the cost function

                  Now we don't necessarily have training data (unless we generate it by using some other method). However, what we do have is the variational principle which allows us to obtain the ground state wave function by minimizing the expectation value of the energy of a trial wavefunction (corresponding to the untrained NQS). Similarly to the traditional variational Monte Carlo method then, it is the local energy we wish to minimize. The gradient to use for the stochastic gradient descent procedure is $$ \begin{align} @@ -1447,7 +599,7 @@ $$











                  -

                  Running the codes

                  +

                  Running the codes

                  @@ -1466,7 +618,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat











                  -

                  Energy as function of iterations, \( N=2 \) electrons

                  +

                  Energy as function of iterations, \( N=2 \) electrons

                  @@ -1477,7 +629,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat











                  -

                  Energy as function of iterations, \( N=6 \) electrons

                  +

                  Energy as function of iterations, \( N=6 \) electrons

                  @@ -1488,7 +640,7 @@ The Broyden-Fletcher-Goldfarb-Shanno algorithm was used to perform the minimizat











                  -

                  Conclusions and where do we stand

                  +

                  Conclusions and where do we stand

                  diff --git a/doc/pub/BM/html/BM.html b/doc/pub/BM/html/BM.html index 42ad4c72f..5524078ad 100644 --- a/doc/pub/BM/html/BM.html +++ b/doc/pub/BM/html/BM.html @@ -6,9 +6,9 @@ Automatically generated HTML file from DocOnce source - + -Solving quantum mechanical problems with Machine Learning and Boltzmann machines +Machine Learning and Boltzmann machines