Monte Carlo Sampling and Estimation of Phase Transition Point#

OpenJij performs simulated annealing (SA). If the temperature is constant, spins can be sampled from the canonical ensemble of that temperature. The canonical ensemble is:

p({σ})=exp(βE({σ}))Z, Z={σ}exp(βE({σ})) p(\{\sigma\}) = \frac{\exp(-\beta E(\{\sigma\}))}{Z}, \ Z = \sum_{\{\sigma\}}\exp(-\beta E(\{\sigma\}))

In the following, we use a fully-connected ferromagnetic Ising model:

E({σ})=JNi<jσiσj (J<0) E(\{\sigma\}) = \frac{J}{N} \sum_{i<j} \sigma_i \sigma_j \ (J<0)

We normalize the Hamiltonian by dividing the energy by the system size NN. Also, we set J=1J = -1.

# import library
import openjij as oj
import numpy as np
import matplotlib.pyplot as plt

# set sampler
sampler = oj.SASampler(num_reads=100)

# define the fully-connected problem
def fully_connected(n):
    h, J = {}, {}
    for i in range(n-1):
        for j in range(i+1, n):
            J[i, j] = -1/n
    return h, J

# set h and J
h, J = fully_connected(n=500)

Let us calculate the magnetization at each temperature.

m=1Niσi m = \frac{1}{N} \sum_i \sigma_i

If the magnetization is closer to 1, the state is in the ordered state (ferromagnetic). Conversely, if it is closer to 0, the state is in the disordered state (paramagnetic). When the temperature is kept constant, it can be seen that the magnetization approaches 0 at a temperature of around 1.0. This is because the phase transition theoretically occurs in the infinite size limit.

# create a list of temperatures
temp_list = np.linspace(0.01, 2, 30)

# calculate the magnetization and its dispersion
mag, mag_std = [], []
for temp in temp_list:
    beta = 1.0/temp
    schedule = [[beta, 100]]
    response = sampler.sample_ising(h, J, schedule=schedule)
    mag_list = [np.abs(np.mean(state)) for state in response.states]
    mag_std.append(np.std(mag_list))
    mag.append(np.mean(mag_list))
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[2], line 9
      7 beta = 1.0/temp
      8 schedule = [[beta, 100]]
----> 9 response = sampler.sample_ising(h, J, schedule=schedule)
     10 mag_list = [np.abs(np.mean(state)) for state in response.states]
     11 mag_std.append(np.std(mag_list))

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/openjij/sampler/sampler.py:222, in BaseSampler.sample_ising(self, h, J, **parameters)
    218 sparse_option = parameters.pop("sparse", True)
    219 bqm = oj.model.model.BinaryQuadraticModel.from_ising(
    220     h, J, sparse=sparse_option
    221 )
--> 222 return self.sample(bqm, sparse=sparse_option, **parameters)

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/openjij/sampler/sa_sampler.py:321, in SASampler.sample(self, bqm, beta_min, beta_max, num_sweeps, num_reads, schedule, initial_state, updater, sparse, reinitialize_state, seed)
    317 sa_system = self._make_system[_updater_name](
    318     _generate_init_state(), ising_graph
    319 )
    320 # ------------------------------------------- choose updater
--> 321 response = self._cxxjij_sampling(
    322     model, _generate_init_state, algorithm, sa_system, reinitialize_state, seed
    323 )
    325 response.info["schedule"] = self.schedule_info
    327 return response

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/openjij/sampler/sampler.py:145, in BaseSampler._cxxjij_sampling(self, model, init_generator, algorithm, system, reinitialize_state, seed, offset)
    140             system_info["system"].append(_sys_info)
    142 # --------------- define sampling execution function
    143 
    144 # Execute sampling function
--> 145 sampling_time = exec_sampling()
    147 # construct response instance
    148 response = oj.sampler.response.Response.from_samples(
    149     states, model.vartype, energies, info=system_info
    150 )

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/openjij/sampler/sampler.py:40, in measure_time.<locals>.wrapper(*args, **kargs)
     37 def wrapper(*args, **kargs):
     38     start_time = time.perf_counter()
---> 40     func(*args, **kargs)
     42     end_time = time.perf_counter()
     43     execution_time = end_time - start_time

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/openjij/sampler/sampler.py:125, in BaseSampler._cxxjij_sampling.<locals>.exec_sampling()
    122     system.reset_spins(init_generator())
    123 # Run sampling algorithm
    124 # and measure execution time
--> 125 _exec_time = measure_time(sampling_algorithm)(system)
    126 execution_time.append(_exec_time)
    128 # get Ising result (state and system information)
    129 # ex. _sys_info save trotterized quantum state.

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/openjij/sampler/sampler.py:40, in measure_time.<locals>.wrapper(*args, **kargs)
     37 def wrapper(*args, **kargs):
     38     start_time = time.perf_counter()
---> 40     func(*args, **kargs)
     42     end_time = time.perf_counter()
     43     execution_time = end_time - start_time

File /opt/hostedtoolcache/Python/3.9.16/x64/lib/python3.9/site-packages/openjij/sampler/sampler.py:99, in BaseSampler._cxxjij_sampling.<locals>.sampling_algorithm(system)
     98 def sampling_algorithm(system):
---> 99     return algorithm(system, self._params["schedule"])

KeyboardInterrupt: 
# plot
plt.errorbar(temp_list, mag, yerr=mag_std)
plt.plot(temp_list, mag)
plt.xlabel('temperature', fontsize=15)
plt.ylabel(r'$|m|$', fontsize=15)
plt.show()
../../../_images/MonteCarloSampling_5_0.png

Phase transition is the phenomenon that changes the value significantly at a specific temperature. In the present model, it is theoretically proven that a phase transition occurs at a temperature of 1.0 in the thermodynamic limit.
However, in many cases, it is difficult to study the properties of the phase transition analytically. For this reason, many methods use Monte Carlo simulations to study the properties of the phase transitions numerically.

Binder Cumulant#

Let us assume that we do not know the phase transition temperature, and let us try to find it as accurately as possible from numerical calculations. The above figure shows that the magnetization is approaching zero as the temperature increases. However, it needs to be made clear which temperature is the phase transition point. This is because the phase transition phenomenon theoretically occurs in a system of infinite size, but the simulation can only handle a finite size, resulting in an error with the theory. This is called the finite-size effect.
It may seem impossible to numerically analyze a system of infinite size. However, in the field of numerical statistical mechanics, methods have been developed to obtain information on the infinite size limit from finite system sizes. One of these methods uses the Binder cumulant. Binder cumulant is:

U412(3m4m22)U_4 \equiv \frac{1}{2}\left( 3- \frac{\langle m^4\rangle}{\langle m^2\rangle^2} \right)
# calculate U_4
def u_4(states):
    m = np.array([np.mean(state) for state in states])
    return 0.5 * (3-np.mean(m**4)/(np.mean(m**2)**2))

We can show that this quantity approaches 1 for ferromagnetic state, where the magnetization approaches 1, and 0 for paramagnetic state, where the magnetization approaches 0. Furthermore, it is known that the value is independent of system size at the phase transition point. Therefore, numerical experiments can be performed for several system sizes, and the phase transition point is the point where U4U_4 for several system sizes intersect. Please check statistical mechanics textbooks for more information. Let us perform the calculation with OpenJij.

# Define the size of the system as a list
n_list = [40, 80, 120, 160]
# Define temperature by list
temp_list = np.linspace(0.5, 1.5, 30)

# Set sampler
sampler = oj.SASampler(num_reads=300)

u4_list_n = []
for n in n_list:
    # Create instance
    h, J = fully_connected(n)
    u4_temp = []
    for temp in temp_list:
        beta = 1.0/temp
        schedule = [[beta, 100 if temp < 0.9 else 300]]
        response = sampler.sample_ising(h, J, 
                                        schedule=schedule, reinitialize_state=False,
                                        num_reads=100 if temp < 0.9 else 1000
                                       )
        u4_temp.append(u_4(response.states))
    u4_list_n.append(u4_temp)
# plot
for n,u4_beta in zip(n_list,u4_list_n):
    plt.plot(temp_list, np.array(u4_beta), label='n={}'.format(n))

plt.legend()
plt.ylabel('$U_4$', fontsize=15)
plt.xlabel('temperature', fontsize=15)
plt.show()
../../../_images/MonteCarloSampling_11_0.png

The statistics are not enough, so the data is a little scattered. However, it can be seen that the data for the four system sizes intersect at a single point around the temperature of 1, which is approximately the phase transition point.
Estimation of the phase transition point by Binder cumulant is a common method used in numerical analysis.

In academic research, it is necessary to obtain sufficient statistics, and of course, error evaluation (calculation of error bars) must be done more precisely. Since the present calculations are only an overview, the exact error evaluations are omitted.

Summary#

We have shown how to perform Monte Carlo sampling using OpenJij. We applied it to an estimation of phase transitions point. OpenJij can be used for various applications depending on your ideas.