Implemented Bose-Hubbard hamiltonian

Created by: kchoo1118

Could you share the script you were using?

Created by: VolodyaCO

@kchoo1118 I have the same issue. Do you know when the abstract machine (with pytorch) and the abstract sampler will be released?

Here is my code:

import netket as nk
import numpy as np
import matplotlib.pyplot as plt
from qutip import destroy, create

def BH_Saito(machine, J=1, U=2, L=11, V=1, n_bosons=9, n_max=5, n_hidden=5):
    
    assert L%2==1
    
    V_list = V * (np.arange(L) - (L-1)//2)**2
    edges = []
    for i in range(L-1):
        edges.append([i, i+1])
        
    g = nk.graph.CustomGraph(edges)
    
    # hilbert space for each site
    hi = nk.hilbert.Boson(graph=g, n_max=n_max, n_bosons=n_bosons)
    
    #boson operators
    a = destroy(n_max + 1).full()
    adag = create(n_max + 1).full()
    num = np.matmul(adag, a)
    
    # fill operators and sites
    operators = []
    sites = []
    
    # Local term
    for i in range(L):
        op_loc = V_list[i] * num + 0.5*U * np.matmul(num, num - np.eye(n_max + 1))
        operators.append(op_loc.tolist())
        sites.append([i])
    
    # Hopping term
    for i in range(L - 1):
        op_loc = -J * np.kron(a, adag)
        operators.append(op_loc.tolist())
        sites.append([i, i+1])
        op_loc = -J * np.kron(a, adag)
        operators.append(op_loc.tolist())
        sites.append([i+1, i])
    
    op = nk.operator.LocalOperator(hi, operators=operators, acting_on=sites)
    
    # Define the fully-connected FFNN    
    if machine == "ffnn":
        layers = (nk.layer.FullyConnected(input_size=L,output_size=n_hidden,use_bias=True),
                  nk.layer.Tanh(input_size=n_hidden),
                  nk.layer.FullyConnected(input_size=n_hidden, output_size=2, use_bias=True),
                  nk.layer.SumOutput(input_size=2))
        for layer in layers:
            layer.init_random_parameters(seed=1234, sigma=0.1)

        ffnn = nk.machine.FFNN(hi, layers)
    
    elif machine == "rbm":
        ffnn = nk.machine.RbmSpinSymm(hi, alpha=n_hidden)
    
    elif machine == "jastrow":
        ffnn = nk.machine.JastrowSymm(hilbert=hi)
        
    elif machine == "bigrbm":
        ffnn = nk.machine.RbmMultiVal(hilbert=hi, alpha=n_hidden)
    
    ffnn.init_random_parameters(seed=0, sigma=0.01)
    
    print('Parameters of wavefunction:', ffnn.n_par)
    
    # Define the sampler that exchanges particles 
    sa = nk.sampler.MetropolisHamiltonian(machine=ffnn, hamiltonian=op)
    #sa = nk.sampler.MetropolisHop(machine=ffnn, d_max=n_bosons)
    
    # Define the optimizer
    opt = nk.optimizer.Sgd(learning_rate=0.03, l2_reg = 1e-5, decay_factor=0.999)
    #opt = nk.optimizer.AdaMax(alpha=0.01)
    #opt = nk.optimizer.AdaDelta()
    
    
    gs = nk.variational.Vmc(hamiltonian=op,
                        sampler=sa,
                        optimizer=opt,
                        n_samples=1000,
                        discarded_samples=200,
                        use_iterative=True,
                        method='Gd')
    
    # Define observables for mean number of bosons in each site:
    for i in range(L):
        op_num_i = [num.tolist()]
        sites_i = [[i]]
        n_i =  nk.operator.LocalOperator(hi, op_num_i, sites_i)
        gs.add_observable(n_i, "n_{0}".format(i))
        
    gs.run(output_prefix='test', n_iter=2000)
    if n_bosons <= 5 and L <= 5:
        res = nk.exact.lanczos_ed(op, first_n=1, compute_eigenvectors=True)
        print("Exact ground state energy = {0:.3f}".format(res.eigenvalues[0]))
        E0 = res.eigenvalues[0]
    else:
        E0 = None
    return ffnn, gs, E0

n_bosons = 4
n_sites = 5
m, vmc, E0 = BH_Saito(machine="bigrbm", L=n_sites, n_max=n_bosons + 1, n_bosons=n_bosons, n_hidden=2)


# Load the data from the .log file
import json

data = json.load(open("test.log"))

# Extract the relevant information
iters = []
energy = []
num = []
for iteration in data["Output"]:
    n = []
    iters.append(iteration["Iteration"])
    energy.append(iteration["Energy"]["Mean"])
    for i in range(n_sites):
        n.append(iteration["n_{0}".format(i)]["Mean"])
    num.append(n)

fig, ax1 = plt.subplots()
ax1.plot(iters, energy, 'k-', label='Energy')
ax1.set_ylabel('Energy')
ax1.set_xlabel('Iteration')
if E0:
    ax1.axhline(float(E0), c='r')
plt.legend()
plt.show()

plt.plot(num[-1], 'ro')
plt.title('Mean occupation of each site')
plt.show()

It is seen that the true ground state is never reached for any of the machines. Also, Saito's work cannot be reproduced with n_bosons=4 and n_sites=11, since the ground energy is never reached (here I do not recommend to use the bigrbm, because the computation time grows a lot. @Mikelesc you can take a look at my code. I reproduced Saito's work with numpy and python correctly, but it is incredibly slow.

Created by: gcarleo

I think it's important to disentangle two things here:

The Hamiltonian, i.e. is the hamiltonian correctly defined?
The variational state, i.e. is a given machine enough to describe efficiently the ground state of that specific hamiltonian

I would first debug 1. and see if the hamiltonian does what we expect, on a small system. To do that, the simplest thing is to call the method to_sparse of the hamiltonian and then use scipy to diagonalize the resulting sparse matrix. If the energies are what you expect for a small system, then we can look into 2 !

Created by: gcarleo

@VolodyaCO

Do you know when the abstract machine (with pytorch) and the abstract sampler will be released?

We are very close, the final release will be in any case very similar to what we have now in the branch called v2.1. If you are interested you can already give it a try. The main thing is that the pytorch module is expected to have two outputs, one for ReLog(Psi) and one for ImLog(Psi)

Created by: VolodyaCO

@gcarleo what worries me is that the full forward neural network should be enough for the problem. I understand it is a FFNN with complex weights. However, it also fails at Saito's problem. Instead, Saito does what one would do with the abstract pytorch machine, which is to define a real NN with two outputs to take into account the phase of the wavefunction.

Created by: gcarleo

@VolodyaCO this strongly depends on what representation is used for the input variables as well... I mean it's an open research problem to understand if it's best to work with binarized inputs over an extended space or just use integers on each physical site

Created by: VolodyaCO

@gcarleo I think that for bosons, the number representation is the best one because the numbers are naturally ordered, whereas in the binarised representation the order needs to be learnt by the machine, adding an extra difficulty. I have tried to install netket 2.1 with python setup.py install, but I'm getting a CMake error:

-- Looking for cheev_ - not found CMake Error at /home/vladimir/.local/share/virtualenvs/Lieb-Liniger-rXmnzkPc/lib/python3.7/site-packages/cmake/data/share/cmake-3.15/Modules/FindLAPACK.cmake:417 (message): A required library with LAPACK API not found. Please specify library location. Call Stack (most recent call first): CMakeLists.txt:238 (find_package)

Created by: gcarleo

You can clearly read the error from there :

A required library with LAPACK API not found. Please specify library

You should install a lapack library. This won't be required in the final version but you need to have it for the moment

Created by: VolodyaCO

@gcarleo I installed them with

sudo apt-get install liblapack-doc liblapack-dev

I leave this if someone else has trouble in the future with this installation. I am running the torch example, which so far has not thrown any errors. @Mikelesc I will introduce this Torch machine to reproduce as closely as possible Saito's paper with netket (although the random sampling with bosons is quite difficult with netket because it relies on spin exchange most of the time).

Created by: VolodyaCO

@Mikelesc I installed netket v2.1 and created the pytorch model to resemble the same neural network used by Saito. The results are bad. Hopefully the input is not binarised with Torch machines (@gcarleo can you confirm this?). If it is not binarised, then I think that the issue is in the sampling. The proposed states have to be randomly generated and accepted/rejected with Metropolis algorithm, and I am not sure that MetropolisHamiltonian sampler allows a "complete" exploration of the configuration space.

Created by: gcarleo

@VolodyaCO no, the input is not binarized. At this stage I'd say this is more of a research related issue rather than netket itself, but here is a couple of thoughts that are relevant also software-wise :

The MetropolisHamiltonian sampler, as detailed in the documentation, uses the hamiltonian matrix elements to generate transitions in the Markov Chain. This means that all symmetries of the hamiltonian are preserved by this sampler
In the case of the Bose Hubbard model, this means that you do a sampling in a subspace at fixed number of particles. In turn, this means that you have to be careful when defining the hilbert object, in providing the correct number of particles you want to find the ground state for
When you say "results are bad", do you mean in comparison to Saito's paper for the very same architecture? This is worth investigating, because it means maybe that the neural network model you are using is not exactly the same?
When comparing to the exact results, one should again take into account the issue of particle number of conservation. Specifically, (and this is indeed an issue with our documentation), our exact diagonalization routine does not preserve the symmetries that are included in the Hilbert object. In the case of the BH model it means that it will find the ground state with the lowest energy, spanning across all possible particle number sectors, that might not be what you are looking for?
I suspect that doing the optimization in the canonical ensemble is much easier than in the grand-canonical and this will give you better energies
I would center the inputs, maybe adding some renormalization/centering layer in a way that the input to the machine is not the bare occupation number but rather n_i - < n >

Created by: VolodyaCO

@gcarleo

I have made the vmc to measure the number of particles in every step, and indeed the number of particles is conserved, so the Hilbert space and the hamiltonian seem to be well defined in my code.
I believe that the machine is the same one as the one defined by Saito. There is an input layer which receives the occupation number at each lattice site, then there is a fully connected hidden layer where a linear superposition of the previous neurons is taken (i.e. the activation function is the identity), and then there is a 2 neuron output layer, where the activation function is tanh. I did this with the torch machine as follows:

model = torch.nn.Sequential(
    torch.nn.Linear(input_size, n_hidden),
    torch.nn.Linear(n_hidden, 2),
    torch.nn.Tanh(),
)

So essentially you are saying that the Lanczos method does not return the ground state for a fixed number of particles?

Thanks for your time

Created by: gcarleo

Indeed ED in netket does not return the energy in the canonical ensemble, so this is important to be considered when comparing to the VMC results!

Created by: VolodyaCO

@gcarleo @Mikelesc I will perform the case with 9 bosons and 11 sites as in Saito's work to check if the energy converges to the same result. If everything goes alright, I think this issue should be closed.

mentioned in issue #330 (closed)

Created by: VolodyaCO

As it can be seen in this image (https://imgur.com/nNCvsBC) the convergence is very poor. RBM goes to ~12, but Saito reports a convergence to the actual ground state (with 9 bosons and 11 sites) with an energy of 10 (in the units of the paper).

Created by: kchoo1118

I think the model should be

model = torch.nn.Sequential( torch.nn.Linear(input_size, n_hidden), torch.nn.Tanh(), torch.nn.Linear(n_hidden, 2), )

Created by: VolodyaCO

I think you can see the results in this image (https://imgur.com/gallery/spnGuDz), and it looks like the results of Saito are practically reached. And it is quite fast, only 30min in my home computer. @Mikelesc

closed

mentioned in merge request !385 (merged)

Implemented Bose-Hubbard hamiltonian

Designs

Child items ...

Activity