Replica exchange Monte Carlo exchange

exchange explores the parameter space by using the replica exchange Monte Carlo (RXMC) method.

Preparation

mpi4py should be installed if MPI parallelization is used.

$ python3 -m pip install mpi4py

Input parameters

This has two subsections algorithm.param and algorithm.exchange .

[algorithm.param]

This defines a space to be explored. When mesh_path key is defined or use_grid key is set to true, the discrete space is used. Otherwise, continuous space is used.

• Continuous space

  • initial_list

    Format: List of float. Length should be equal to dimension.

    Description: Initial value of parameters. If not defined, these will be initialized randomly.

  • min_list

    Format: List of float. Length should be equal to dimension.

    Description: The minimum value that each parameter can take.

  • max_list

    Format: List of float. Length should be equal to dimension.

    Description: The maximum value that each parameter can take.

  • step_list

    Format: List of float. Length should be equal to dimension.

    Description: Step width (standard deviation of the Gaussian distribution) used for Monte Carlo updates.

  • pbc_list

    Format: List of boolean (default: false for each parameter).

    Description: Whether to use periodic boundary conditions (PBC) for generating local candidate points for each parameter. The length should be equal to dimension.

  • About state proposals

    • In Monte Carlo updates, candidate points are proposed according to a Gaussian distribution centered at the current state. If a candidate parameter d falls outside the range [min_list[d], max_list[d]):

      • If pbc_list[d] = true, it will be wrapped into the range [min_list[d], max_list[d]) using periodic boundary conditions.

      • If pbc_list[d] = false, the proposal will be rejected.

• Discrete space

  • mesh_path

    Format: string

    Description: Path to the mesh definition file. See Reference file below for the format.

  • comments

    Format: String (default: “#”)

    Description: Character(s) that indicate the beginning of a comment line when reading the mesh definition file.

  • delimiter

    Format: String (default: whitespace)

    Description: Column delimiter of the mesh definition file. Specify "," to read a CSV file.

  • skiprows

    Format: Integer (default: 0)

    Description: Number of lines to skip at the beginning of the mesh definition file. Use it to skip header lines.

  • neighborlist_path

    Format: string

    Description: Path to the neighborhood-list file. See Reference file below for the format. When it is omitted, the neighborhood list will be generated automatically assuming that the mesh points within radius are neighbors.

  • radius

    Format: Float.

    Description: the points within the distance specified by radius are considered as neighbors. It is mandatory when neighborlist_path is omitted or use_grid is set to true.

  • use_grid

    Format: Boolean.

    Description: When it is set to true, a uniform mesh will be generated using the parameters min_list, max_list, and num_list.

  • min_list

    Format: List of float. Length should be equal to dimension.

    Description: Lower end of each parameter of the mesh to be generated.

  • max_list

    Format: List of float. Length should be equal to dimension.

    Description: Upper end of each parameter of the mesh to be generated.

  • num_list

    Format: List of integer. Length should be equal to dimension.

    Description: The number of mesh points along each parameter of the mesh to be generated.

[algorithm.exchange]

• numsteps

  Format: Integer

  Description: The number of Monte Carlo steps.

• numsteps_exchange

  Format: Integer

  Description: The number of Monte Carlo steps between replica exchange trials.

• numsteps_thermalization

  Format: Integer

  Description: The number of Monte Carlo steps to discard for thermalization. The default is 10% of numsteps.

• Tmin

  Format: Float

  Description: The minimum value of the “temperature” (\(T\)).

• Tmax

  Format: Float

  Description: The maximum value of the “temperature” (\(T\)).

• bmin

  Format: Float

  Description: The minimum value of the “inverse temperature” (\(\beta = 1/T\)). Specify either the temperature range (Tmin, Tmax) or the inverse-temperature range (bmin, bmax), but not both pairs.

• bmax

  Format: Float

  Description: The maximum value of the “inverse temperature” (\(\beta = 1/T\)). Specify either the temperature range (Tmin, Tmax) or the inverse-temperature range (bmin, bmax), but not both pairs.

• Tlogspace

  Format: Boolean (default: true)

  Description: Whether to assign “temperature” to replicas equally spaced in the logarithmic space or not.

• nreplica_per_proc

  Format: Integer (default: 1)

  Description: The number of replicas in an MPI process. The total number of replicas is given by “number of MPI processes × nreplica_per_proc”. In the replica exchange Monte Carlo method, one temperature point is assigned to each replica, and therefore the number of temperature points generated between Tmin and Tmax equals the total number of replicas. To use more temperature points, increase the number of MPI processes or nreplica_per_proc.

• separate_T

  Format: Boolean (default: true)

  Description: Whether to write log files of Monte Carlo steps separately for each temperature (result_T*.txt). This option is ignored when export_combined_files is true.

• export_combined_files

  Format: Boolean (default: false)

  Description: Whether to write the contents of trial.txt and result.txt into a single combined file combined.txt instead of separate per-process files. Use the odatse_extract_combined tool to extract the individual files from the combined file (see odatse_extract_combined).

Note

For developers: Setting the environment variable ODATSE_USE_MPI_BUFFERED=1 switches the MPI communication used for collecting replica data from object-based (gather) to buffer-based (Gather). It is normally not needed, but is provided as a performance-tuning option for large-scale parallel runs.

Reference file

Mesh definition file

Define the grid space to be explored in this file. The first column is the index of the mesh, and the second and subsequent columns are the values of variables. Note that the index of the mesh will be ignored for this “algorithm”.

Below, a sample file is shown.

1 6.000000 6.000000
2 6.000000 5.750000
3 6.000000 5.500000
4 6.000000 5.250000
5 6.000000 5.000000
6 6.000000 4.750000
7 6.000000 4.500000
8 6.000000 4.250000
9 6.000000 4.000000
...

Neighborhood-list file

Before searching a discrete space by the Markov chain Monte Carlo method, the “neighborhood” of each point \(i\) must be defined, i.e. the points to which a walker can move from \(i\). A neighborhood-list file defines these neighborhoods. In this file, the index of an initial point \(i\) is specified by the first column, and the indices of final points \(j\) are specified by the second and successive columns.

An utility tool, odatse_neighborlist is available for generating a neighborhood-list file from a mesh file. For details, please see Related Tools.

0 1 2 3
1 0 2 3 4
2 0 1 3 4 5
3 0 1 2 4 5 6 7
4 1 2 3 5 6 7 8
5 2 3 4 7 8 9
...

Output files

In the following files, the unit that tracks the state of each replica is called a walker (nreplica_per_proc walkers per process).

RANK/trial.txt

This file stores the parameter candidates proposed (trialed) in the Monte Carlo updates and the corresponding values of the objective function for each replica. The first column is the index of the MC step. The second column is the index of the walker in the process. The third column is the temperature of the replica. The fourth column is the value of the solver. The remaining columns are the coordinates.

Example:

# step walker T fx z1 z2
0 0 0.004999999999999999 0.07830821484593968 3.682008067401509 3.9502750191292586
1 0 0.004999999999999999 0.0758494287185766 2.811346329442423 3.691101784194861
2 0 0.004999999999999999 0.08566823949124412 3.606664760390988 3.2093903670436497
3 0 0.004999999999999999 0.06273922648753057 4.330900869594549 4.311333132184154

RANK/result.txt

This file stores the sequence of accepted parameters after the acceptance/rejection step and the corresponding values of the objective function for each replica (when a trial is rejected, the previous state is kept). This has the same format as trial.txt.

# step walker T fx z1 z2
0 0 0.004999999999999999 0.07830821484593968 3.682008067401509 3.9502750191292586
1 0 0.004999999999999999 0.07830821484593968 3.682008067401509 3.9502750191292586
2 0 0.004999999999999999 0.07830821484593968 3.682008067401509 3.9502750191292586
3 0 0.004999999999999999 0.06273922648753057 4.330900869594549 4.311333132184154

best_result.txt

The optimal value of the solver and the corresponding parameter among all the samples.

nprocs = 4
rank = 2
step = 65
walker = 0
fx = 0.008233957976993406
z1 = 4.221129370933539
z2 = 5.139591716517661

result_T#.txt

This file stores samples for each temperature (# is replaced with the index of the temperature). The first column is the index of the MC step. The second column is the index of the walker. The third column is the value of the solver. The remaining columns are the coordinates.

# T = 1.0
0 15 28.70157662892569 3.3139009347685118 -4.20946994566609
1 15 28.70157662892569 3.3139009347685118 -4.20946994566609
2 15 28.70157662892569 3.3139009347685118 -4.20946994566609
3 15 28.98676409223712 3.7442621319489637 -3.868754990884034

fx.txt

This file stores statistical metrics over all the replicas for each temperature. The first column is the inverse temperature. The second and third columns are the expectation value and the standard error of the solver’s output (\(f(x)\)), respectively. The fourth column is the number of replicas. The fifth column is the logarithm of the ratio between the normalization factors (partition functions)

\[\log\frac{Z}{Z_0} = \log\int \mathrm{d}x e^{-\beta f(x)} - \log\int \mathrm{d}x e^{-\beta_0 f(x)},\]

where \(\beta_0\) is the minimum value of \(\beta\) used in the calculation. The sixth column is the acceptance ratio of MC updates.

# $1: 1/T
# $2: mean of f(x)
# $3: standard error of f(x)
# $4: number of replicas [Not used for exchange MC]
# $5: log(Z/Z0)
# $6: acceptance ratio
0.01 71.1757607526179 0.7074887893238192 80 0.0 0.906
0.011236548001387516 71.64077624599895 0.6352374919714137 80 -0.08478181677306895 0.9014444444444445
0.012626001098748564 65.8095785630431 0.6095581760655332 80 -0.18096066826346802 0.8936666666666667
0.014187266741165962 60.51980595731712 0.5843304334378878 80 -0.2798580927722558 0.8766666666666667
  ...

Restart

The execution mode is specified by the run_mode parameter to the constructor. The operation of each mode is described as follows. The parameter values correspond to --init, --resume, and --cont options of odatse command, respectively.

• "initial" (default)

  The program is started from the initialized state. If the checkpointing is enabled, the intermediate states will be stored on the following occasions:

  1. the specified number of steps has been performed, or the specified period of time has passed.

  2. at the end of the execution.

• "resume"

  The program execution is resumed from the latest checkpoint. The conditions such as the number of MPI processes should be kept the same.

• "continue"

  The program execution is continued from the previous run. The value of numsteps should be increased. The step counter is taken over.

  For example: in the first run, the calculation is carried out for 1000 steps with numsteps = 1000. In the next run, the calculation is continued with numsteps = 2000, where the calculations from 1001st step to 2000th step are carried out.

Algorithm

Markov chain Monte Carlo

The Markov chain Monte Carlo (MCMC) sampling explores the parameter space by moving walkers \(\vec{x}\) stochastically according to the weight function \(W(\vec{x})\). For the weight function, the Boltzmann factor \(W(\vec{x}) = e^{-f(\vec{x})/T}\) is generally adopted, where \(T>0\) is the “temperature.” Unfortunately, in many cases it is impossible to sample walkers directly according to \(W\). Instead, the MCMC method moves walkers by small steps and generates a time series \(\{\vec{x}_t\}\) such that the distribution of the walkers obeys \(W\) . Let us denote the transition probability from \(\vec{x}\) to \(\vec{x}'\) by \(p(\vec{x}' | \vec{x})\). When \(p\) is determined by the following condition (“the balance condition”)

\[W(\vec{x}') = \sum_{\vec{x}} p(\vec{x}' | \vec{x}) W(\vec{x}),\]

the distribution of the generated time series \(\{\vec{x}_t\}\) will converge to \(W(\vec{x})\) [1]. Practically, the stronger condition (“the detailed balance condition”)

\[p(\vec{x} | \vec{x}') W(\vec{x}') = W(\vec{x})p(\vec{x}' | \vec{x})\]

is usually imposed. The detailed balance condition reduces to the balance condition when summed over \(\vec{x}\).

ODAT-SE adopts the Metropolis-Hastings (MH) method for solving the detailed balance condition. The MH method splits the transition process into the proposal process and the acceptance process.

1. Generate a candidate \(\vec{x}\) with the proposal probability \(P(\vec{x} | \vec{x}_t)\).

   • As \(P\), use a simple distribution such as a normal distribution centered at \(\vec{x}_t\).

2. Accept the candidate \(\vec{x}\) with the acceptance probability \(Q(\vec{x} | \vec{x}_t)\).

   • If accepted, let \(\vec{x}_{t+1}\) be \(\vec{x}\).

   • Otherwise, let \(\vec{x}_{t+1}\) be \(\vec{x}_t\).

The whole transition probability is the product of these two, \(p(\vec{x} | \vec{x_t}) = P(\vec{x} | \vec{x}_t) Q(\vec{x} | \vec{x}_t)\). The acceptance probability \(Q(\vec{x} | \vec{x}_t)\) is defined as

\[Q(\vec{x} | \vec{x}_t) = \min\left[1, \frac{W(\vec{x})P(\vec{x}_t | \vec{x}) }{W(\vec{x}_t) P(\vec{x} | \vec{x}_t)} \right].\]

It is easy to verify that the detailed balance condition is satisfied by substituting it into the detailed balance condition equation.

When adopting the Boltzmann factor for the weight and a symmetric distribution \(P(\vec{x} | \vec{x}_t) = P(\vec{x}_t | \vec{x})\) for the proposal probability, the acceptance probability \(Q\) takes the following simple form:

\[Q(\vec{x} | \vec{x}_t) = \min\left[1, \frac{W(\vec{x})}{W(\vec{x}_t)} \right] = \min\left[1, \exp\left(-\frac{f(\vec{x}) - f(\vec{x}_t)}{T}\right) \right].\]

By defining \(\Delta f = f(\vec{x}) - f(\vec{x}_t)\) and using the fact \(Q=1\) for \(\Delta f \le 0\), the procedure of MCMC with the MH algorithm is the following:

1. Choose a candidate from near the current position and calculate \(f\) and \(\Delta f\).

2. If \(\Delta f \le 0\), that is, the walker is descending, accept it.

3. Otherwise, accept it with the probability \(Q=e^{-\Delta f/T}\).

4. Repeat 1-3.

The solution is given as the point giving the minimum value of \(f(\vec{x})\). The third process of the above procedure ensures that walkers can climb over hills of height \(\Delta f \sim T\), so that the MCMC sampling can escape from local minima.

Replica exchange Monte Carlo

The “temperature” \(T\) is one of the most important hyperparameters in the MCMC sampling. The MCMC sampling can climb over hills of height \(T\) but cannot easily escape from valleys deeper than \(T\). This is why the temperature should be increased in order to avoid getting stuck in local minima. On the other hand, since walkers cannot resolve valleys smaller than \(T\), the precision of the obtained result \(\min f(\vec{x})\) becomes about \(T\), and it is necessary to decrease the temperature in order to achieve a more precise result. This dilemma means that the temperature must be tuned carefully.

One way to overcome this problem is to update the temperature as well. For example, simulated annealing decreases the temperature as the iterations proceed. Another algorithm, simulated tempering, treats temperature as another parameter to be sampled, not a fixed hyper parameter, and updates the temperature after some iterations according to the (detailed) balance condition. Simulated tempering studies the details of a valley by cooling and escapes from a valley by heating. Replica exchange Monte Carlo (RXMC), also known as parallel tempering, is a parallelized version of simulated tempering. In this algorithm, several copies of the system at different temperatures, called replicas, are simulated in parallel. Then, with some interval of steps, each replica exchanges temperature with another one according to the (detailed) balance condition. As the simulated tempering does, RXMC can observe the details of a valley and escape from it by cooling and heating. Moreover, because each temperature is assigned to just one replica, the temperature distribution will not be biased. Using more replicas narrows the temperature interval, and increases the acceptance ratio of the temperature exchange. This is why this algorithm is well suited to massively parallel calculations.

It is recommended that users perform minsearch optimization starting from the result of exchange, because the RXMC result has uncertainty due to temperature.

footnote

[1]

To be precise, aperiodicity and ergodicity are necessary for convergence.