Population Annealing Monte Carlo `pamc`#

pamc explores the parameter space by using the Population Annealing Monte Carlo (PAMC) method.

Preparation#

mpi4py should be installed for the MPI parallelization.

$ python3 -m pip install mpi4py

Input parameters#

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

`[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.pamc]`#

numsteps

Format: Integer

Description: The number of Monte Carlo steps.
numsteps_annealing

Format: Integer

Description: The number of Monte Carlo steps between successive temperature reductions.
Tnum

Format: Integer

Description: The number of “temperature” points.
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 (population size) is given by “number of MPI processes × nreplica_per_proc”.
resampling_interval

Format: Integer (default: 1)

Description: The number of annealing steps between resamplings of the replicas.
fix_num_replicas

Format: Boolean (default: true)

Description: Whether to fix the number of replicas or not on resampling.
separate_T

Format: Boolean (default: true)

Description: Whether to write log files of Monte Carlo steps separately for each temperature. 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, result.txt, and weight.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).
anneal_from_beta0

Format: Boolean (default: false)

Description: When set to true and bmin>0, the algorithm first performs annealing and resampling from an initial random sample at \(\beta=0\) (infinite temperature) up to the smallest specified inverse temperature (bmin or \(1/T_{\max}\)), before starting the main calculation. This ensures that, even if the computation does not start from \(\beta=0\), the reference value of \(\log Z/Z_0\) can still be considered as being taken at \(\beta=0\).

About the number of steps#

Specify just two of numsteps, numsteps_annealing, and Tnum. The value of the remaining one will be determined automatically. They are approximately related by numsteps = numsteps_annealing × Tnum (when the division leaves a remainder, the extra steps are assigned to the higher-temperature points).

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.

6.000000 6.000000
6.000000 5.750000
6.000000 5.500000
6.000000 5.250000
6.000000 5.000000
6.000000 4.750000
6.000000 4.500000
6.000000 4.250000
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.

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

Output files#

`RANK/trial_T#.txt`#

This file stores the suggested parameters and the corresponding value returned from the solver for each temperature point (specified by #). The first column (step) is the index of the MC step. The second column (walker) is the index of the walker in the process. The third column (beta) is the inverse temperature of the replica (when the temperature range is given via Tmin/Tmax, the temperature T is written instead). The fourth column (fx) is the value of the solver. The fifth through (4+dimension)-th columns are the coordinates. The last two columns (weight and ancestor) are the Neal-Jarzynski weight and the grand-ancestor of the replica.

Example:

# step walker beta fx x1 weight ancestor
0 0.0 73.82799488298886 8.592321856342956 1.0 0
1 0.0 13.487174782058675 -3.672488908364282 1.0 1
2 0.0 39.96292704464803 -6.321623766458111 1.0 2
3 0.0 34.913851603463 -5.908794428939206 1.0 3
4 0.0 1.834671825646121 1.354500581633733 1.0 4
5 0.0 3.65151610695736 1.910894059585031 1.0 5
...

`RANK/trial.txt`#

This file combines all the trial_T#.txt files into one.

`RANK/result_T#.txt`#

This file stores the sampled parameters and the corresponding value returned from the solver for each replica and each temperature. This has the same format as trial_T#.txt.

# step walker beta fx x1 weight ancestor
0 0.0 73.82799488298886 8.592321856342956 1.0 0
1 0.0 13.487174782058675 -3.672488908364282 1.0 1
2 0.0 39.96292704464803 -6.321623766458111 1.0 2
3 0.0 34.913851603463 -5.908794428939206 1.0 3
4 0.0 1.834671825646121 1.354500581633733 1.0 4
5 0.0 3.65151610695736 1.910894059585031 1.0 5
...

`RANK/result.txt`#

This file combines all the result_T#.txt files into one.

`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

`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 (when anneal_from_beta0 = true, the reference is \(\beta_0 = 0\)). 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
# $5: log(Z/Z0)
# $6: acceptance ratio
0.0 33.36426034198166 3.0193077565358273 100 0.0 0.9804
0.1 4.518006242920819 0.9535301415484388 100 -1.2134775491597027 0.9058
0.2 1.5919146358616842 0.2770369776964151 100 -1.538611313376179 0.9004
...

`RANK/weight.txt`#

This file stores the Neal-Jarzynski weight of each replica at each temperature. The columns are, in order, the temperature index (Tindex), the inverse temperature (beta), the walker index (walker), the grand-ancestor id (idnum), the value of the objective function (fx), the logarithm of the weight (log_weight), and the coordinates.

Example:

# Tindex beta walker idnum fx log_weight x1
0 0.0 0 0 73.82799488298886 0.0 8.592321856342956
0 0.0 1 1 13.487174782058675 0.0 -3.672488908364282
...

`pr.txt`#

This file stores the participation ratio of the weights for each temperature, which measures the effective number of replicas. The first column is the temperature index (Tindex), the second column is the inverse temperature (\(1/T\)), and the third column is the participation ratio.

# $1: Tindex
# $2: 1/T
# $3: participation ratio
0 0.0 100.0
1 0.1 87.23456789012345
...

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. when the calculation at a temperature point finishes and 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 sequence of the temperature points should be specified so that it is continuous from that of the previous run.

If the temperature was lowered from Tmax= \(T^{(1)}\) to Tmin= \(T^{(2)}\) in the previous run, the next values should be taken as Tmax= \(T^{(2)}\) and Tmin= \(T^{(3)}\). In the new calculation, the temperature points are taken from \(T^{(2)}\) to \(T^{(3)}\) divided by Tnum, namely, \(T_0 = T^{(2)}\), \(T_1\),…, \(T_{\text{Tnum}-1}=T^{(3)}\). (Tnum can be different from the previous run.)

Algorithm#

Goal#

When the weight of the configuration \(x\) under some parameter \(\beta_i\) is given as \(W_i(x)\) (e.g., the Boltzmann factor \(W_i(x) = \exp[-\beta_i f(x)]\) ), the expectation value of \(A\) is defined as

\[\langle A\rangle_i = \frac{\int \mathrm{d}xA(x)W_i(x)}{\int \mathrm{d}x W_i(x)} = \frac{1}{Z_i}\int \mathrm{d}xA(x)W_i(x) = \int \mathrm{d}xA(x)\tilde{W}_i(x),\]

where \(Z_i = \int \mathrm{d} x W_i(x)\) is the normalization factor (partition function) and \(\tilde{W}_i(x) = W_i(x)/Z_i\) is the probability of \(x\).

Our goal is to numerically calculate the expectation value for each \(\beta_i\) and the (ratios of the) normalization factors.

Annealed Importance Sampling (AIS) [1]#

First, we introduce a series of configurations \(\{x_i\}\) obeying the following joint probability

\[\tilde{W}(x_0, x_1, \dots, x_n) = \tilde{W}_n(x_n) \tilde{p}_n(x_n, x_{n-1}) \tilde{p}_{n-1}(x_{n-1}, x_{n-2}) \cdots \tilde{p}_1(x_1, x_0),\]

with

\[\tilde{p}_i(x_i, x_{i-1}) = p_i(x_{i-1}, x_i) \frac{\tilde{W}_i(x_{i-1})}{\tilde{W}_i(x_i)},\]

where \(p_i(x, x')\) is a transition probability from \(x\) to \(x'\) under \(\beta_i\) holding the balance condition,

\[\int \mathrm{d}x \tilde{W}_i(x) p_i(x, x') = \tilde{W}_i(x').\]

It turns out that \(\tilde{W}_n(x_n)\) is the marginal distribution of \(\tilde{W}(x_0, x_1, \dots, x_n)\), that is,

\[\tilde{W}_n(x_n) = \int \prod_{i=0}^{n-1} \mathrm{d} x_i \tilde{W}(x_0, x_1, \dots, x_n),\]

from

\[\int \mathrm{d} x_{i-1} \tilde{p}_i(x_i, x_{i-1}) = \int \mathrm{d} x_{i-1} \tilde{W}_i(x_{i-1}) p_i(x_{i-1}, x_i) / \tilde{W}_i(x_i) = 1.\]

Consequently, \(\langle A \rangle_n\) is represented by using the extended configuration \(\{x_i\}\) as

\[\begin{split}\begin{split} \langle A \rangle_n &\equiv \int \mathrm{d} x_n A(x_n) \tilde{W}_n(x_n) \\ &= \int \prod_i \mathrm{d} x_i A(x_n) \tilde{W}(x_0, x_1, \dots, x_n). \end{split}\end{split}\]

Unfortunately, it is difficult to generate directly a series of configurations \(\{x_i\}\) following the distribution \(\tilde{W}(x_0, x_1, \dots, x_n)\). Then, instead of \(\tilde{W}(x_0, x_1, \dots, x_n)\), we consider \(\{x_i\}\) obeying the joint distribution

\[\tilde{g}(x_0, x_1, \dots, x_n) = \tilde{W}_0(x_0) p_1(x_0, x_1) p_2(x_1, x_2) \dots p_n(x_{n-1}, x_n),\]

by using the following scheme:

Generate \(x_0\) from the initial distribution \(\tilde{W}_0(x)\)
Generate \(x_{i+1}\) from \(x_i\) through \(p_{i+1}(x_i, x_{i+1})\)

By using the reweighting method (or importance sampling method), \(\langle A \rangle_n\) is rewritten as

\[\begin{split}\begin{split} \langle A \rangle_n &= \int \prod_i \mathrm{d} x_i A(x_n) \tilde{W}(x_0, x_1, \dots, x_n) \\ &= \int \prod_i \mathrm{d} x_i A(x_n) \frac{\tilde{W}(x_0, x_1, \dots, x_n)}{\tilde{g}(x_0, x_1, \dots, x_n)} \tilde{g}(x_0, x_1, \dots, x_n) \\ &= \left\langle A\tilde{W}\big/\tilde{g} \right\rangle_{g, n} \end{split}.\end{split}\]

Because the ratio between \(\tilde{W}\) and \(\tilde{g}\) is

\[\begin{split}\begin{split} \frac{\tilde{W}(x_0, \dots, x_n)}{\tilde{g}(x_0, \dots, x_n)} &= \frac{\tilde{W}_n(x_n)}{\tilde{W}_0(x_0)} \prod_{i=1}^n \frac{\tilde{p}_i(x_i, x_{i-1})}{p_i(x_{i-1}, x_i)} \\ &= \frac{\tilde{W}_n(x_n)}{\tilde{W}_0(x_0)} \prod_{i=1}^n \frac{\tilde{W}_i(x_{i-1})}{\tilde{W}_i(x_i)} \\ &= \frac{Z_0}{Z_n} \frac{W_n(x_n)}{W_0(x_0)} \prod_{i=1}^n \frac{W_i(x_{i-1})}{W_i(x_i)} \\ &= \frac{Z_0}{Z_n} \prod_{i=0}^{n-1} \frac{W_{i+1}(x_{i})}{W_i(x_i)} \\ &\equiv \frac{Z_0}{Z_n} w_n(x_0, x_1, \dots, x_n), \end{split}\end{split}\]

the form of the expectation value will be

\[\langle A \rangle_n = \left\langle A\tilde{W}\big/\tilde{g} \right\rangle_{g, n} = \frac{Z_0}{Z_n} \langle Aw_n \rangle_{g,n}.\]

Finally, the ratio between the normalization factors \(Z_n/Z_0\) can be evaluated as

\[\frac{Z_n}{Z_0} = \langle w_n \rangle_{g,n},\]

and therefore the expectation value of \(A\) can be evaluated as a weighted arithmetic mean:

\[\langle A \rangle_n = \frac{\langle Aw_n \rangle_{g,n}}{\langle w_n \rangle_{g,n}}.\]

This weight \(w_n\) is called the Neal-Jarzynski weight.

Population annealing (PA) [2]#

Although the AIS method can estimate the expectation values of \(A\) for each parameter \(\beta\) in the form of a weighted arithmetic mean, the variance of the weights \(w\) is generally large, so the accuracy of the result deteriorates. In order to overcome this problem, the population annealing Monte Carlo (PAMC) method periodically resamples all the replicas according to the probability \(p^{(k)} = w^{(k)} / \sum_k w^{(k)}\) and resets all the weights to unity.

The following pseudo code describes the scheme of PAMC:

for k in range(K):
    w[0, k] = 1.0
    x[0, k] = draw_from(β[0])
for i in range(1, N):
    for k in range(K):
        w[i, k] = w[i-1, k] * ( W(x[i-1,k], β[i]) / W(x[i-1,k], β[i-1]) )
    x_prev = x[i-1, :]
    if i % interval == 0:
        x_prev = resample(x_prev, w[i, :])
        w[i, :] = 1.0
    for k in range(K):
        x[i, k] = transfer(x_prev[k], β[i])
    a[i] = sum(A(x[i,:]) * w[i,:]) / sum(w[i,:])

There are two resampling methods: one with a fixed number of replicas [2] and one without [3].

References#

[1] R. M. Neal, Statistics and Computing 11, 125-139 (2001).

[2] K. Hukushima and Y. Iba, AIP Conf. Proc. 690, 200 (2003).

[3] J. Machta, Phys. Rev. E 82, 026704 (2010).

Population Annealing Monte Carlo pamc#

Preparation#

Input parameters#

[algorithm.param]#

[algorithm.pamc]#