# CNN Encoding of Acoustic Parameters for Prominence Detection

###### Abstract

Expressive reading, considered the defining attribute of oral reading fluency,
comprises the prosodic realization of phrasing and prominence. In the context
of evaluating oral reading, it helps to establish the speaker’s comprehension
of the text. We consider a labeled dataset of children’s reading recordings
for the speaker-independent detection of prominent words using acoustic-
prosodic and lexico-syntactic features. A previous well-tuned random forest
ensemble predictor is replaced by an RNN sequence classifier to exploit
potential context dependency across the longer utterance. Further, deep
learning is applied to obtain word-level features from low-level acoustic
contours of fundamental frequency, intensity and spectral shape in an end-to-
end fashion. Performance comparisons are presented across the different
feature types and across different feature learning architectures for
prominent word prediction to draw insights wherever possible.

Index Terms: word prominence, children’s speech, oral reading evaluation

##  1 Introduction

The prosodic structure of speech carries important information in terms of the
syntax and the meaning, both of which are critical to a listener’s ease of
comprehension of the spoken message [1, 2, 3]. Phrase boundaries embed
sentence syntax through word grouping while prominence or emphasis on specific
words signals new information or highlights a contrast. In this paper, we
investigate methods for the automatic detection of prominent words in the oral
reading of middle-school children in the context of second language learning.
The appropriate use of prosody is indicative of the reader’s comprehension of
the text and hence constitutes a critical component of oral reading evaluation
systems [4, 5, 6]. In beginning readers, suprasegmental skills take longer to
develop compared to word decoding ability, with phrasing coming earlier than
the effective use of prominence.

Prominence is perceived by a listener when a word stands out of its local
context in one or more of the suprasegmental attributes such as duration, F0,
intensity and spectral shape [7]. The local context itself refers to the
phones and syllables within the word as well as a neighborhood of up to
several words. Prosody perception is influenced not only by the low-level
acoustic cues but also top-down expectations from lexico-syntactic information
[8, 9, 10]. The precise combination and relative importance of the cues
depends on the speaker, language and speaking style as also on the listener.
Various aggregates of the sampled acoustic parameters across the word segment
including mean and variance, contour shape descriptors, and differences in
these quantities across neighboring words comprise word-level prosodic
features. These features are then used to train a conventional supervised
classifier for the automatic detection, possibly in combination with lexico-
syntactic information [11, 12, 13, 14]. In our own recent work, we applied
systematic feature selection within and across the distinct suprasegmental
attributes in a random forest ensemble predictor to derive a compact set of
interpretable features for speaker-independent boundary and prominence
detection on a children’s oral reading dataset [15]. With the search space for
such ‘hand-crafted’ features being very large, however, the process can miss
potentially important features. Further, the pre-selected context windows used
in such analyses make it difficult to exploit the long and variable time scale
of prosodic relationships across an utterance in any comprehensive manner. The
potential for deep learning solutions has therefore been recognized for some
time but incorporated successfully in the prominence detection task only more
recently, as briefly reviewed next.

Rosenberg et al. [16] used a large number of acoustic-prosodic features and
aggregates at word level derived from their previous work [11, 17] (which gave
rise to the AuToBI tool) in a BiRNN classifier where the word sequence context
was learned over that explicitly provided in the feature vector. They observed
a small improvement ($<$ 1% absolute) in boundary and pitch accent detection
over a baseline conditional random forest classifier. Wu et al. [18] also used
similar aggregated acoustic features in an LSTM to find an improvement over
the use of an SVM classifier. Lin et al. [19] used a hierarchical BLSTM
network to aggregate features across phone, syllable and word to model
contextual information at multiple granularities in the joint detection of
boundaries and prominence.

In a departure from pre-computed word-level features, a recent work [20, 21]
uses CNN on frame-based acoustic parameters (energy, F0, loudness, voicing
probability, zero crossing rate and harmonic-to-noise ratio) together with a
context window of two neighboring words to optimally learn the high-level
aggregation features. The max-pooled CNN feature maps are directly classified
with a softmax layer. With word position indicators provided in the input
segment, they report an improvement of 1-3% points absolute over Rosenberg
[11] on lexical stress and phrase boundary detection on the BURNC corpus, with
speaker-independent scenarios being more challenging. Zhang et al. [22] also
use acoustic contours and MFCC features over 10 s segments as inputs to a CNN
with fixed narrow kernel width of 3 frames, with syllable and word position
indicators marked at the frame level. The CNN outputs go to a BLSTM classifier
to obtain emphasis at frame level.

Given the significance of both local acoustic features and longer, more
global, contexts spanning several words and possibly different sentences
across the utterance in the perception of prominence, it is reasonable to
consider architectures combining low-level feature aggregation with sequence
models. The CNN learned feature representation of Stehwien et al. [20] was
extended recently by Nielsen et al. [23] using full utterances as input and
adding an LSTM layer to incorporate greater context. An improvement in
accuracy of 1% was noted for pitch accent detection on the BURNC corpus, with
a further similar increase when they added text features by concatenating
various dimension Glove word embeddings with the CNN embedding at the LSTM
input.

Inspired by the above reviewed works, we investigate specific enhancements to
our previously proposed random forest based prominence prediction for
children’s oral reading evaluation, with its highly optimized acoustic-
prosodic word level features [15, 24]. The same dataset of children’s read
stories in English is used with its 42,138 words across 800 utterances by 35
speakers, recorded at 16 kHz sampling and manually transcribed at word level.
The selected speakers have reasonable word decoding ability in English (as
second language) but widely varying levels of prosodic skill. The individual
utterances comprise between 50-70 words, each word labeled for the
presence/absence of prominence by 7 naive listeners using RPT [25]. With a
binary prominence decision based on 3 or more votes, we obtain a reasonable
figure of 24% of the total words labeled prominent. We begin with replacing
the random forest predictor with an RNN sequence model that can, in principle,
capture implicit context dependence from across the utterance. Next, lexico-
syntactic features based on the linguistic correlates of prominence are
included, with new features related to the canonical structure of the text.
Finally, motivated by the demonstrated potential of convolutional neural
networks to learn discriminative patterns and thus replace any feature
engineering, we investigate CNN architectures for our task in the end-to-end
learning of prosodic word-level features from low-level acoustic parameters.

##  2 RNN-based prominence scoring

Using utterances segmented via forced alignment with the transcript, Sabu and
Rao [15] obtained a compact and highly optimized set of word-level features by
applying random forest model based feature selection on the children’s speech
dataset considered here. A very large set of acoustic-prosodic features
defined across the different suprasegmental attributes of pitch, duration,
intensity and spectral balance were considered, with multiple ways of defining
temporal context in a $\pm$2 word neighborhood, to obtain a reduced set of 34
features. A Pearson correlation of 0.69 was obtained between the random forest
regression output and the degree of prominence based on proportion of rater
votes, and detection F-score of 0.63, in the speaker-independent prominent
word detection task. In this section, we report work on using the same set of
acoustic-prosodic features with an RNN model with its input given by the
variable length sequence of words across an input utterance. We consider
further the inclusion of lexico-syntactic information in the input sequence.

###  2.1 Architecture

We tried different RNN architectures: GRU [26] and LSTM [27] in both
unidirectional and bidirectional configurations. The number of layers were
varied in the range {1, 2, 3} while the number of hidden units were picked
from the range {48, 96, 128, 256, 512}. At each time-step, a feature vector
corresponding to a single word is fed as input to the model. A final feed-
forward layer linearly transforms the RNN output at each time step to a
scalar, which is passed through a sigmoid layer to get the degree of
prominence prediction.

###  2.2 Adding lexical features

The high predictive power of lexico-syntactic information, such as part-of-
speech tags, has been noted for several corpora before [28, 29]. Baumann [10]
investigated a wide range of non-prosodic factors for prominence prediction in
German speech with random forest based feature selection. Lexical information
included 12 PoS tags: NN (noun), NP (proper noun), JJ (adjective), RB
(adverb), VB (verb), AU (auxiliary verb), MD (modal verb), PR (pronoun), IN
(preposition), CC (conjunction), RP (particle) and DT (article). Assuming that
this fine break-up may not be suited to our dataset of learners who do not
necessarily have high levels of text comprehension, we also consider more
coarse groupings of different dimensions. A one-hot encoding is used for the
tags per word with the PoS itself determined from English Grammar rules [30]
since automatic parsers showed poor performance. We also include the number of
phones and syllables per word which indirectly relates to word frequency [10].

We further propose the canonical information structure (i.e. expected
prominence and phrase boundaries) as additional useful features linked to top-
down cues. We determine this by applying syntax and givenness rules [3].
Motivated by the observation that expected prosodic events depend on reading
speed, the events are labeled mandatory, optional and forbidden. With no known
NLP methods for the automatic extraction of this information from arbitrary
text, we use the model reading of the story to validate our labels. With
reading miscues, albeit few, being a part of our dataset, the PoS and
information structure tags are based on the target word after achieving the
automatic alignment of the transcript with the text. Other lexical features
such as the number of phones and number of syllables in the word are estimated
for the uttered word rather and inserted words are marked with a uniform PoS
tag.

##  3 Learning word-level features with CNN

The word-level features of our previous work were computed from word (and sub-
word) aligned contours corresponding to the time-varying acoustic parameters
of F0, intensity and spectral shape, computed at 10 ms intervals [15].
Utterance based z-score normalization is applied on F0 and intensity
parameters. We wish to investigate CNN-based automatic learning of word-level
features from the same low-level acoustic contours. Given the previously
observed speaker-dependence of the relative importances of the different
prosodic attributes, we investigate a 3-channel CNN architecture where
attribute-wise embeddings are computed with their own best filters and
concatenated for the final representation [22]. The contour groups are F0 (4
contours), intensity (4 contours) and spectral shape including HNR and
spectral band energies (7 contours) and each feature group is input to
separate CNN filter bank as shown in Figure 1. The output encodings from the
filter banks are finally concatenated and fed to the sequence classifier
(GRU/LSTM) network together with other word-level features as considered in
our experiments.

Perceived prominence depends on the underlying temporal dynamics while CNN
computations are time-invariant. To match the information available to the
hand-crafted feature extraction, the CNN input is a segment with the selected
word context, possibly with word and subword positions within the segment
provided explicitly. We explore a range of context choices as also different
types of position information, starting with the context of $\pm$ 1 word given
the importance of immediately neighboring words in the realization of emphasis
on a word [15, 21]. We also consider CNN filter banks with a range of kernel
width choices motivated by sub-word units [31]. The 1D convolution output of
each of the CNN filters is max-pooled across time to get a scalar value per
filter per filter bank. Each filter bank has N filters, each with k kernel
sizes, resulting in a kN-dimensional feature encoding, for each channel,
corresponding to a word.

Figure 1: Prominence degree prediction architecture.

##  4 Experiments and results

###  4.1 Train-test splits and baseline

The random forest predictor based system reported in [15] serves as a baseline
and we adopt the identical training and testing methodology here. The complete
dataset of 42,138 words is split into three equal folds with no speaker
overlap for 3-fold cross-validation based testing. The hyperparameters are
tuned with 4-fold CV on the train split. The model is then trained on the
entire train set and results are reported on the corresponding unseen test
set. We report the mean and standard deviation across the three test folds.
The results for the prominence degree prediction are reported in terms of
Pearson correlation between the predictor output and the degree of prominence
from the RPT rater votes. Also, the prominence detection F-score values are
reported considering prominence present when a word receives 3 or more votes
as discussed in Section 1.

###  4.2 RNN training and performance

We present the entire input utterance in the form of a sequence of word-level
features to an RNN classifier. Various feature combinations are investigated
to determine the individual and combined contributions of prosodic and lexical
features. We test various RNN architectures as presented in Section 2. For
training, we used the AdamW [32] optimizer, an improved version of popular
Adam [33] optimizer with a weight decay mechanism that helps with faster
training and more generalized models. We used a learning rate of 0.003 and a
batch size of 500. Dropout with probability 0.2 is added to each RNN layer
except the last layer. The Mean Squared Error between the scaled target score
(between 0 and 1) and the predicted score is used as the loss function to be
minimized during training.

Table 1: Performance of various models with set of 34 acoustic features. (* indicates sd $<$ 0.01) Model | # layers | # units | Correlation | F-score
---|---|---|---|---
RFC | - | - | 0.69* | 0.63*
GRU | 2 | 96 | 0.68 | 0.63
LSTM | 2 | 256 | 0.69 | 0.63
BGRU | 2 | 96 | 0.70 | 0.64
BLSTM | 2 | 256 | 0.71* | 0.64*

Table 1 shows the performance with the 34 acoustic-prosodic features (‘A34’)
used in the RNN architectures compared with the same features in the baseline
random forest predictor. We note that there is an improvement from the
utterance-long context available to the sequence models, especially in the
case of the bidirectional models. We employ the BGRU going ahead due to it
ease of training and known suitability for lower dataset sizes.

Next, considering separately the lexical features, Figure 2 shows performance
for various reduced sets starting from the full set of 21 that includes PoS
tags and phone/syllable counts (together termed ‘L’), and information
structure labels (termed ‘I’). From a maximum correlation of 0.72, we note a
drop in performance as the 6 information structure and 2 word length features
are removed to get to 13 features corresponding to the PoS alone. A further
drop is recorded when the PoS tags are grouped in different ways to get to the
final feature set of only content/function word distinctions. We note that
lexical features alone show a predictive power similar to prosodic features
alone (and not much higher as in some previous work [18, 29]), explained by
the lower proficiency speakers of our dataset. All the same, the reduction of
PoS tags clearly hurts performance. Table 2 summarizes the achieved
performance gains as we augment the acoustic-prosodic feature set with lexical
(all but information structure) and information structure features with each
of the latter two clearly adding value.

Figure 2: Correlation performance variation with different sets of lexical features. Table 2: Performance with addition of lexical and information structure features. (* indicates sd $<$ 0.01) Features | Correlation | F-score
---|---|---
A34 | 0.70 | 0.64
A34 + L | 0.75* | 0.67*
A34 + L + I | 0.79* | 0.69*

###  4.3 CNN training and performance

While our multi-channel CNN framework is similar to that of Zhang et al. [22],
we expand the search for architecture choices by considering the use of
multiple kernel widths in each channel to capture the distinct time scales of
acoustic variation. We start from the 4 kernels with widths [5, 11, 25, 51]
similar to that of the sentence parsing CNN architecture of Trang et al. [31],
which roughly cover sub-phone, phone, syllable and word, and possibly some
context. Given the fixed narrow kernel width of 3 frames used in [20, 22], we
add this to our candidates for testing. From the different combinations
presented in Table 3, we observe that the syllable and word width kernel sizes
(25, 51) helps the performance while including other widths does not change
it. For these instances, we fixed the number of filters of each kernel size to
8. The input context to CNN is also fixed to $\pm$ 1 word with position
encoding used to indicate current word based on [21]. We find that phone and
sub-phone width kernels do not help, and even degrade the performance in some
cases. We finalize the use of two kernel widths 25 frames and 51 frames,
corresponding to syllable and word widths, for each CNN filter bank. Next, to
find the optimal number of kernels for each width, we varied the number of CNN
filters in the range {4, 8, 16, 20}. We found that 8 filters give the best
performance and adopted this for all the further experiments.

Table 3: Performance for different sets of kernel width choices with the corresponding lit. reference indicated. Kernel widths | Correlation | F-score
---|---|---
3,5,11,25,51 | 0.67 | 0.62
5,11,25,51 [31] | 0.67 | 0.62
25,51 | 0.67* | 0.62*
11,25,51 | 0.67* | 0.62
11 [23] | 0.65 | 0.61*

Next, we consider different choices for position encoding. Stehwien et al.
[20] found it useful to indicate the frames corresponding to the word as
distinct from the context frames from the neighboring words within the input
segment. However, it has been observed that the preceding and following words
influence perceived prominence in different ways [15]. Therefore, we decided
to change the position encoding to distinguish all the three words. To test
whether a special encoding to indicate inter-word pauses can help, we applied
5-bit one-hot encoding across the segment of 3 words and 2 pauses. The results
are shown in Table 4. Taking position information to intra-word level, we
explore one-hot encoding to indicate syllable number (up to a maximum of 7
syllables) in the word as indicated in last row of the table. As can be seen,
the pause encoding improves the performance thereby indicating that it is
important indeed to differentiate between next and previous word as well as
pauses. Syllable marking helps too, and further experiments are reported using
the 5-bit encoding along with the 7-bit syllable position.

Table 4: Performance with various positional encoding approaches. Input to CNN has $\pm$ 1 word context. Positional encoding | Correlation | F-score
---|---|---
1-bit [21] | 0.67* | 0.62*
3-bit (word order) | 0.68 | 0.62
5-bit (word/pause order) | 0.68 | 0.62
5-bit + syllable position | 0.68 | 0.63

We also experimented with choosing fixed duration input segments, centered at
the word, as an alternative to the variable (3 word) segments considered so
far to avoid the zero-padding at the CNN input. We find that the duration that
works best corresponds to the average 3-word segment and that there is no
performance gain with fixing input segment length.

To capture the contribution of multi-channel processing of the attribute-wise
contour, we compare it with single-channel processing of the combined contours
in Table 5 with all else kept unchanged. We note a drop in performance.
Feature pooling of multi-channel outputs also reduced performance compared to
simple feature concatenation suggesting the importance of retaining all
attribute variations in the input to the RNN model.

Table 5: Performance with different CNN architectures Architecture | Correlation | F-score
---|---|---
Multichannel (concatenate) | 0.69 | 0.63
Multichannel (pooling) | 0.66 | 0.62
Single-channel | 0.67* | 0.62

###  4.4 Overall performance

Given our overall goal of investigating the automatic learning of features for
a prosodic prominence detection task on a challenging dataset of children’s
read speech, we explored the cascade of CNN and RNN with word and sentence
level inputs respectively. Our baseline was a random forest ensemble predictor
with hand-crafted acoustic-prosodic features optimized with to exploit the
best of acoustic parameter aggregation, normalization and local context for
prominence detection of a word within a long utterance. Using a bidirectional
GRU with the sequence of acoustic-prosodic features (‘A34’) helped improve
F-score by 1% over the baseline. As we see in row 2 of Table 6, the CNN-
learned features in the same setting fall slightly short.

Adding features assumed to be harder to learn automatically, but important for
word prominence, such as the actual and speech-rate normalized word, syllable
and pause durations, we concatenate a 12-dimensional duration feature vector
(‘DP-12’) with the CNN-learned features. We further select 10 contour shape
features (‘A10’, a subset of A34) that require fitting the temporal variation
of the acoustic parameters. The added features are seen to bring distinct
additional value. Finally, the inclusion of text features (L, I) boosts
performance for both feature sets as we see in the last 2 rows of Table 6 with
higher performance for A34.

Our outcomes are overall consistent with those of the few previous works that
have compared automatically learned and hand-crafted features on the same
dataset and task. As is well known, automatic feature learning is limited by
dataset sizes and future work must examine the use of larger, possibly
unlabeled, data and architectures that learn additional useful information
such as speaker identity or employ attention for better context learning.

Table 6: Performance of CNN encoding concatenated with different word-level features as RNN input. (* indicates sd $<$ 0.01) Features | Correlation | F-score
---|---|---
A34 | 0.70 | 0.64
CNN | 0.69 | 0.63
CNN + D-P12 + A10 | 0.71 | 0.64
CNN + D-P12 + A10 + L + I | 0.77* | 0.68
A34 + L + I | 0.79* | 0.69*

## References

  * [1] J. K. Bock and J. R. Mazzella, “Intonational marking of given and new information: Some consequences for comprehension,” _Memory & Cognition_, vol. 11, no. 1, pp. 64–76, 1983.
  * [2] L. van Maastricht, T. Zee, E. Krahmer, and M. Swerts, “L1 perceptions of L2 prosody: The interplay between intonation, rhythm, and speech rate and their contribution to accentedness and comprehensibility,” in _Proceedings of INTERSPEECH_ , Stockholm, Sweden, 2017, pp. 364–368.
  * [3] J. M. Levis and A. O. Silpachai, “Prominence and information structure in pronunciation teaching materials,” in _Proceedings of the Pronunciation in Second Language Learning and Teaching conference_ , Ames, IA, USA, 2017, pp. 216–229.
  * [4] S. E. Sinambela, “Prosody as a tool for assessing reading fluency of adult ESL students,” _Advances in Language and Literary Studies_ , vol. 8, no. 6, pp. 83–87, 2017.
  * [5] D. Paige, W. Rupley, G. Smith, W. N. T. Rasinski, and T. Magpuri-Lavell, “Is prosodic reading a strategy for comprehension?” _Journal for Educational Research_ , vol. 141, no. 4, pp. 245–275, 2017.
  * [6] H. K. Keskin, G. Ari, and M. Bastug, “Role of prosodic reading in listening comprehension,” _International Journal of Education & Literacy Studies_, vol. 7, no. 1, pp. 59–65, 2019.
  * [7] M. Breen, E. Fedorenko, M. Wagner, and E. Gibson, “Acoustic correlates of information structure,” _Language and Cognitive Processes_ , vol. 25, no. 7/8/9, pp. 1044–1098, 2010.
  * [8] J. Cole, Y. Mo, and M. Hasegawa-Johnson, “Signal-based and expectation-based factors in the perception of prosodic prominence,” _Laboratory Phonology_ , vol. 1, no. 2, pp. 425–452, 2010.
  * [9] P. Wagner, F. Tamburini, and A. Windmann, “Robust tracking for automatic reading tutors,” in _Proceedings of INTERSPEECH_ , Portland, OR, USA, 2012\.
  * [10] S. Baumann and B. Winter, “What makes a word prominent? Predicting untrained German listeners’ perceptual judgments,” _Journal of Phonetics_ , vol. 70, pp. 20–38, 2018.
  * [11] A. Rosenberg, “Automatic detection and classification of prosodic events,” Ph.D. dissertation, Columbia University, 2009.
  * [12] T. Mishra, V. R. Sridhar, and A. Conkie, “Word prominence detection using robust yet simple prosodic features,” in _Proceedings of INTERSPEECH_ , Portland, OR, USA, 2012, pp. 1864–1867.
  * [13] G. Christodoulides and M. Avanzi, “An evaluation of machine learning methods for prominence detection in french,” in _Proceedings of INTERSPEECH_ , Singapore, 2014, pp. 116–119.
  * [14] M. P. Black, D. Bone, Z. I. Skordilis, R. Gupta, W. Xia, P. Papadopoulos, S. N. Chakravarthula, B. Xiao, M. V. Segbroeck, J. Kim, P. G. Georgiou, and S. S. Narayanan, “Automated evaluation of non-native english pronunciation quality: Combining knowledge- and data-driven features at multiple time scales,” in _Proceedings of INTERSPEECH_ , Dresden, Germany, 2015, pp. 493–497.
  * [15] K. Sabu and P. Rao, “Prosodic event detection in children’s read speech,” _Computer Speech and Language_ , vol. 68, pp. 1–19, 2021.
  * [16] A. Rosenberg, R. Fernandez, and B. Ramabhadran, “Modeling phrasing and prominence using deep recurrent learning,” in _Proceedings of INTERSPEECH_ , Dresden, Germany, 2015, pp. 3066–3070.
  * [17] A. Rosenberg, “AuToBI - A tool for automatic ToBI annotation,” in _Proceedings of INTERSPEECH_ , Makuhari, Japan, 2010, pp. 146–149.
  * [18] Y. Wu, H. Li, and S. Li, “Automatic pitch accent detection using long short-term memory neural networks,” in _Proceedings of International Symposium on Signal Processing Systems_ , Beijing, China, 2019.
  * [19] B. Lin, L. Wang, X. Feng, and J. Zhang, “Joint detection of sentence stress and phrase boundary for prosody,” in _Proceedings of INTERSPEECH_ , Shanghai, China, 2020, pp. 4392–4396.
  * [20] S. Stehwien and N. T. Vu, “Prosodic event recognition using convolutional neural networks with context information,” in _Proceedings of INTERSPEECH_ , Stockholm, Sweden, 2017, pp. 2326–2330.
  * [21] S. Stehwien, A. Schweitzer, and N. T. Vu, “Acoustic and temporal representations in convolutional neural network models of prosodic events,” _Speech Communication_ , vol. 125, pp. 128–141, 2020.
  * [22] L. Zhang, F. M. Jia Jia and, S. Zhou, W. Chen, C. Zhang, and R. Li, “Emphasis detection for voice dialogue applications using multi-channel convolutional bidirectional long short-term memory network,” in _Proceedings of INTERSPEECH_ , Hyderabad, India, 2018.
  * [23] E. Nielsen, M. Steedman, and S. Goldwater, “The role of context in neural pitch accent detection in English,” in _Proceedings of Conference on Empirical Methods in Natural Language Processing_ , 2020.
  * [24] K. Sabu and P. Rao, “Detection of prominent words in oral reading by children,” in _Proceedings of Speech Prosody_ , Poznan, Poland, 2018, pp. 314–318.
  * [25] J. Cole, T. Mahrt, and J. Roy, “Crowd-sourcing prosodic annotation,” _Computer Speech and Language_ , vol. 45, pp. 300–325, 2017.
  * [26] J. Chung, Ç. Gülçehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” _CoRR_ , vol. abs/1412.3555, 2014. [Online]. Available: http://arxiv.org/abs/1412.3555
  * [27] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural Computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [28] K. Chen, M. H.-J. A. Cohen, and J. Cole, “A maximum likelihood prosody recognizer,” in _Proceedings of Speech Prosody_ , Nara, Japan, 2004.
  * [29] S. Ananthakrishnan and S. S. Narayanan, “Automatic prosodic event detection using acoustic, lexical, and syntactic evidence,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 16, pp. 216–228, 2008.
  * [30] P. Wren, H. Martin, and N. P. Rao, _High School English Grammar and Composition_. New Delhi, India: S. Chand & Company Ltd., 1936.
  * [31] T. Tran, S. Toshniwal, M. Bansal, K. Gimpel, K. Livescu, and M. Ostendorf, “Parsing speech: A neural approach to integrating lexical and acoustic-prosodic information,” in _Proceedings of NAACL-HLT_ , New Orleans, Louisiana, 2018.
  * [32] S. Gugger and J. Howard, “Adamw and super-convergence is now the fastest way to train neural nets,” https://www.fast.ai/2018/07/02/adam-weight-decay/, 2018, accessed: 2021-04-01.
  * [33] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2017.


# Exchanging replicas with unequal cost, infinitely and permanently

Sander Roet  Daniel T. Zhang  Titus S. van Erp<EMAIL_ADDRESS>Department of Chemistry, Norwegian University of Science and Technology
(NTNU),
N-7491 Trondheim, Norway

###### Abstract

We developed a replica exchange method that is effectively parallelizable even
if the computational cost of the Monte Carlo moves in the parallel replicas
are considerably different, for instance, because the replicas run on
different type of processor units or because of the algorithmic complexity. To
prove detailed-balance, we make a paradigm shift from the common conceptual
viewpoint in which the set of parallel replicas represents a high-dimensional
superstate, to an ensemble based criterion in which the other ensembles
represent an environment that might or might not participate in the Monte
Carlo move. In addition, based on a recent algorithm for computing permanents,
we effectively increase the exchange rate to infinite without the steep
factorial scaling as function of the number of replicas. We illustrate the
effectiveness of the replica exchange methodology by combining it with a
quantitative path sampling method, replica exchange transition interface
sampling (RETIS), in which the costs for a Monte Carlo move can vary
enormously as paths in a RETIS algorithm do not have the same length and the
average path lengths tend to vary considerably for the different path
ensembles that run in parallel. This combination, coined $\infty$RETIS, was
tested on three model systems.

Replica Exchange $|$ Path Sampling $|$ infinite swapping $|$ Markov-chain
Monte Carlo

The Markov chain Monte Carlo (MC) method is one of the most important
numerical techniques for computing averages in high-dimensional spaces, like
the configuration space of a many particle system. The approach has
applications in a wide variety of fields ranging from computational physics,
theoretical chemistry, economics, and genetics. The MC algorithm effectively
generates a selective random walk through state space in which the artificial
steps are designed such to ensure that the frequency of visiting any
particular state is proportional to the equilibrium probability of that state.
The Metropolis [1] or the more general Metropolis-Hastings [2] algorithms are
the most common approaches for designing such random steps (MC moves) based on
the detailed-balance principle. That is, the MC moves should be constructed
such that the number of transition from an old state $s^{(o)}$ to a new state
$s^{(n)}$ is exactly balanced by the number of transitions from the new to the
old state: $\rho(s^{(o)})\pi(s^{(o)}\rightarrow
s^{(n)})=\rho(s^{(n)})\pi(s^{(n)}\rightarrow s^{(o)})$ where, $\rho(\cdot)$ is
the state space equilibrium probability density and $\pi(\cdot)$ are the
probabilities to make a transition between the two states given the set of
possible MC moves. Further, the transition is split into a generation and an
acceptance/rejection step such that $\pi(s\rightarrow s^{\prime})=P_{\rm
gen}(s\rightarrow s^{\prime})P_{\rm acc}(s\rightarrow s^{\prime})$. In the
case that the sampled state space is the configuration space of a molecular
system at constant temperature, $P_{\rm gen}$ might relate to moving a
randomly picked particle in a random direction over a small random distance,
and $\rho(s)$ is proportional to the Boltzmann weight $e^{-\beta E(s)}$ with
$\beta=1/k_{B}T$ the inverse temperature and $E(s)$ the state’s energy. The
Metropolis-Hastings algorithm takes a specific solution for the acceptance
probability

$\displaystyle P_{\rm acc}(s^{(o)}\rightarrow s^{(n)})={\rm
min}\left[1,\frac{\rho(s^{(n)})P_{\rm gen}(s^{(n)}\rightarrow
s^{(o)})}{\rho(s^{(o)})P_{\rm gen}(s^{(o)}\rightarrow s^{(n)})}\right]$ (1)

The generation probabilities will cancel in the above expression if they are
symmetric, $P_{\rm gen}(s\rightarrow s^{\prime})=P_{\rm
gen}(s^{\prime}\rightarrow s)$ as in the less generic Metropolis scheme. At
each MC step, the new state is either accepted or rejected based on the
probability above. In case of a rejection, the old state is maintained and
resampled. This scheme obeys detailed-balance and if, in addition, the set of
MC moves are ergodic, equilibrium sampling is guaranteed. When ergodic
sampling, even if mathematically obeyed, is slowed down by a rough (free)
energy landscape, Replica exchange MC becomes useful.

Replica exchange MC (or replica exchange molecular dynamics) is based on the
idea to simulate several copies of the system with different ensemble
definitions [3, 4, 5], most commonly ensembles with increasing temperature
(parallel tempering). By performing “swaps” between adjacent replicas, the
low-temperature replicas gain access to the broader space region that are
explored by the high-temperature replicas. The detailed-balance and
corresponding acceptance-rejection step can be derived by viewing the set of
states in the different ensembles (replicas) as a single high-dimensional
superstate $S=(s_{1},s_{2},\ldots,s_{N})$ representing the system in a set of
$N$ independent ”parallel universes”. The Metropolis scheme applied to the
superstate yields

$\displaystyle P_{\rm acc}(S^{(o)}\rightarrow S^{(n)})={\rm
min}\left[1,\frac{\rho(S^{(n)})}{\rho(S^{(o)})}\right]$ (2)

in which the probability of the superstate equals

$\displaystyle\rho(S)=\rho(s_{1},s_{2},\ldots,s_{N})=\prod_{i=1}^{N}\rho_{i}(s_{i})$
(3)

where $\rho_{i}(\cdot)$ is the specific probability density of ensemble $i$.
For example, the move that attempts to swap the first two states, implying
$S^{(o)}=(s_{1},s_{2},\ldots,s_{N})$ and $S^{(n)}=(s_{2},s_{1},\ldots,s_{N})$,
will be accepted with a probability

$\displaystyle P_{\rm acc}={\rm
min}\left[1,\frac{\rho_{1}(s_{2})\rho_{2}(s_{1})}{\rho_{1}(s_{1})\rho_{2}(s_{2})}\right]$
(4)

In a replica exchange simulation, swapping moves and standard MC or MD steps
are applied alternately. Parallel computing will typically distribute the same
number of processing units per ensemble to carry out the computational
intensive standard moves. The swapping move is cheap, but it requires that the
ensembles involved in the swap have completed their previous move. If the
standard moves in each ensemble require different computing times, then
several processing units have to wait for the slow ones to finish. If the
disbalance per move is relatively constant, the replicas could effectively be
made to progress in cohort by trying to differentiate the number of processing
units per ensemble or the relative frequency of doing replica exchange versus
standard moves per ensemble. However, in several MC methods this disbalance is
not constant, such as with configurational bias MC [6, 7, 8] or path sampling
[9]. The number of elementary steps to grow a polymer in configurational bias
MC obviously depends on the polymer’s length that is being grown, but also
early rejections lead to a broad distribution of the time it takes to complete
a single MC move even in uniform polymer systems. Analogously, the time
required to complete a MC move in path sampling simulations will depend on the
length of the path being created. Other examples of complex Monte Carlo
methods with a fluctuating CPU cost per move are cluster Monte Carlo
algorithms [10] and event-chain Monte Carlo [11, 12].

We will show that the standard acceptance Eqs. 1 and 4 can be applied in a
parallel scheme in which ensembles are updated irregularly in time and the
average frequency of MC moves is different for the ensembles. In addition, we
show that we can apply an infinite swapping [13] scheme between the available
ensembles. For this, we develop a new protocol based on the evaluation of
permanents that circumvents the steep factorial scaling. This last development
is also useful for standard replica exchange.

## Methods

Finite swapping. In the following, we will assume that we have two types of MC
moves. One move that is CPU intensive and can be carried out within a single
ensemble, and replica exchange moves between ensembles which are relatively
cheap to execute. The CPU intensive move will be carried out by a single
worker (one processor unit, one node or a group of nodes) and these workers
perform their task in parallel on the different ensembles. One essential part
of our algorithm is that we have less workers than ensembles such that
whenever the worker is finished and produced a new state for one ensemble,
this state can directly be swapped with the states of any of the available
ensembles (the ones not occupied by a worker). After that, the worker will
randomly switch to another unoccupied ensemble for performing a CPU intensive
move.

In its most basic form, the algorithm consists of the following steps:

  1. 1.

Define $N$ ensembles and let $\rho_{i}(\cdot)$ be the probability distribution
of ensemble $i$. We also define $P_{\rm RE}$ which is the probability for a
replica exchange move.

  2. 2.

Assign $K<N$ ’workers’ (processing units) to $K$ of the $N$ ensembles for
performing a CPU intensive MC move. Each ensemble is at all times occupied by
either 1 or 0 workers. The following steps are identical for all the workers.

  3. 3.

If the worker is finished with its MC move in ensemble $i$, the new state is
accepted or rejected according to Eq. 1 (with $\rho_{i}$ for $\rho$). Ensemble
$i$ is updated with the new state (or by resampling the old state in case of
rejection) and is then considered to be free.

  4. 4.

Take a uniform random number $\nu$ between 0 and 1. If $\nu>P_{\rm RE}$ go to
step 7.

  5. 5.

Among the available ensembles, pick a random pair $(i,j)$.

  6. 6.

Try to swap the states of ensembles $i$ and $j$ using Eq. 4 (with labels $i,j$
instead of $1,2$). Update ensembles $i,j$ with the swapped state or the old
state in case of a rejection. Return to step 4.

  7. 7.

Select one of the free ensembles at random and assign the worker to that
ensemble for performing a new standard move. Go to step 3.

In this algorithm ensembles are not updated in cohort like in standard replica
exchange, but updates occur at irregular intervals. In addition, the different
ensemble conditions can result in systematic differences in the number of
states that are being created over time. To prove that the above scheme
actually samples the correct distributions requires a fundamentally new
conceptual view as the superstate picture is no longer applicable. Despite
that the algorithm uses the same type of Eqs. 1 and 4, as one would use in
standard replica exchange, it does not rely on Eqs. 2 and 3 that are no longer
valid. In the Supplementary Information (SI) we provide a proof from the
individual ensemble’s perspective in which the other ensembles provide an
”environment” ${\mathcal{E}}$ that might, or might not, participate in the
move of the ensemble considered. By doing so, we no longer require that the
number of transitions from old to new, $S^{(o)}\rightarrow S^{(n)}$, is the
same as from new to old, $S^{(n)}\rightarrow S^{(o)}$. Instead, by writing
$S=(s_{1},{\mathcal{E}})$, from ensemble 1’s perspective, we have that the
number of
$(s_{1}^{(o)},{\mathcal{E}}^{(o)})\rightarrow(s_{1}^{(n)},\textrm{${}^{a}\hskip-2.0pt{\mathcal{E}}$}^{(n)})$
transitions should be equal to the number of
$(s_{1}^{(n)},{\mathcal{E}}^{(o)})\rightarrow(s_{1}^{(o)},\textrm{${}^{a}\hskip-2.0pt{\mathcal{E}}$}^{(n)})$
transitions when the standard move is applied where
$\textrm{${}^{a}\hskip-2.0pt{\mathcal{E}}$}^{(n)}$ refers to any new
environment. The SI shows a similar detailed-balance condition for the replica
exchange moves. At step 6 we sample only ensemble $i$ and $j$ or,
alternatively, all free ensembles get a sample update. This would mean
resampling the existing state of those not involved in a swap (”null move”).
This makes the approach more similar to the superstate sampling albeit using
only free ensembles, as described in the SI. The null move does not reduce the
statistical uncertainty, but we mention it here as it makes it easier to
explain the infinite swapping approach. But for the detailed-balance
conditions to be valid it is imperative that occupied ensembles are not
sampled.

An essential aspect of the efficiency of our algorithm is that the number of
workers $K$ is less than the number of ensembles $N$. The case $K=N$ is valid
but would reduce the number of replica exchange moves to zero as only one
ensemble is free at the maximum. Reducing the $K/N$ ratio will generally imply
a higher acceptance in the replica exchange moves as we can expect a higher
number of free ensembles whose distributions have significant overlap. What
gives the optimum number of workers is therefore a non-trivial question that
we will further explore in the Results and Discussion section. However, for
case $K<N$ we can maximize the effect of the replica exchange moves by taking
the $P_{\rm RE}$ parameter as high as possible. In fact, we can simulate the
effect of the limit $P_{\rm RE}\rightarrow 1$ without having to do an infinite
number of replica exchange moves explicitly. This lead to an infinite swapping
[13] version of our algorithm.

Infinite swapping. If in the previously described algorithm we take $P_{\rm
RE}=1-\delta$, we will loop through the steps 4-6 for many iterations ($n_{\rm
it}=\sum_{n=0}^{\infty}n(1-\delta)^{n}\delta=1/\delta$ in the limit
$\delta\rightarrow 0$) before getting to step 7. When $\delta$ vanishes and
$n_{\rm it}$ becomes infinitely large, we expect that all possible swaps will
be executed an infinite number of times. Since the swaps obey detailed balance
between unoccupied ensembles, these will essentially sample the distribution
of Eq. 3 (for the subset $S^{*}$ of unoccupied ensembles). Hence, when the
loop is exited, each possible permutation $\sigma\in S^{*}$ has been sampled
$n_{\rm it}\times\rho(\sigma)/\sum_{\sigma}\rho(\sigma)$ times. By lumping all
the times that the same permutation was sampled and normalizing by division
with $n_{\rm it}$, we simply sample all the possible permutations in one go
using fractional weights that sum up to 1. This is then the only sampling
step, as the single update in step 3 can be skipped due to its negligible
$1/n_{\rm it}$ weight.

The idea of doing an ”infinite number” of swapping moves has been proposed
before [13, 14, 15], but here we give a different flavor to this approach by a
convenient reformulation of the problem into permanents that allows us to beat
the steep factorial scaling reported in earlier works [13]. The permanents
formulation goes as follows. Supposed that after step 3, there are 4 free
ensembles (we name them $e_{1},e_{2},e_{3},e_{4}$) containing 4 states
($s_{1}$, $s_{2}$, $s_{3}$, $s_{4}$). Which state is in which ensemble after
this step is irrelevant. We can now define a weight-matrix $W$:

$W=\bordermatrix{&e_{1}&e_{2}&e_{3}&e_{3}\cr
s_{1}&W_{11}&W_{12}&W_{13}&W_{14}\cr s_{2}&W_{21}&W_{22}&W_{23}&W_{24}\cr
s_{3}&W_{31}&W_{32}&W_{33}&W_{34}\cr s_{4}&W_{41}&W_{42}&W_{43}&W_{44}}$

where $W_{ij}\propto\rho_{j}(s_{i})$. Essential to our approach is the
computation of the permanent of the $W$ matrix, ${\rm perm}(W)$, and that of
the $W{\\{ij\\}}$-matrices in which the row $i$ and column $j$ are removed.

The permanent of a matrix is similar to the determinant, but without
alternating signs. We can, henceforth, write ${\rm
perm}(W)=\sum_{j=1}^{4}W_{1j}{\rm perm}(W{\\{1j\\}})$. As the permanent of the
$1\times 1$ matrix is obviously equal to the single matrix value, the
permanent of arbitrary dimension could in principle be solved recursively
using this relation. Based on the permanents of $W$, we will construct a
probability matrix $P$:

$P=\bordermatrix{&e_{1}&e_{2}&e_{3}&e_{4}\cr
s_{1}&P_{11}&P_{12}&P_{13}&P_{14}\cr s_{2}&P_{21}&P_{22}&P_{23}&P_{24}\cr
s_{3}&P_{31}&P_{32}&P_{33}&P_{34}\cr s_{4}&P_{41}&P_{42}&P_{43}&P_{44}}$

where $P_{ij}$ is the chance to find state $s_{i}$ in ensemble $e_{j}$. As for
each permutation each state is in one ensemble and each ensemble contains one
state, the $P$-matrix is bistochastic: both the columns and the rows sum up to
1. If we consider $S^{*}_{ij}$ the set of permutations in which state $s_{i}$
is in $e_{j}$, we can write $P_{ij}=\sum_{\sigma\in
S^{*}_{ij}}\rho(\sigma)/\sum_{\sigma^{\prime}\in S^{*}}\rho(\sigma^{\prime})$.
We can, however, also use the permanent representation in which

$\displaystyle P_{ij}=\frac{W_{ij}{\rm perm}(W\\{ij\\})}{{\rm perm}(W)}$ (5)

So far we have not won anything as computing the the permanent via the
recursive relation mentioned above has still the factorial scaling. The
Gaussian elimination approach, that allows an order ${\mathcal{O}}(n^{3})$
computation for determinants of $n\times n$ matrices, won’t work for
permanents as only some but not all row- and column-operations have the same
effect to a permanent as to a determinant. One can for instance swap rows and
columns without changing the permanent. Multiplying a row by a nonzero scalar
multiplies the permanent by the same scalar. Hence, this will not affect the
P-matrix based on Eq. 5. Unlike the determinant, adding or subtracting to a
row a scalar multiple of another row, an essential part of the Gaussian
elimination method, does change the permanent. This makes the permanent
computation of a large matrix excessively more expensive than the computation
of a determinant. Yet, recent algorithms based on the
Balasubramanian–Bax–Franklin–Glynn (BBFG) formula [16, 17, 18] scale as
${\mathcal{O}}(2^{n})$. This means that the computation of the full $P$-matrix
scales as ${\mathcal{O}}(2^{n}\times n^{2})$, which seems still steep but is
nevertheless a dramatic improvement compared to factorial scaling.

For our target time of 1 second, for instance, we could only run the algorithm
up to $N=7$ in the factorial approach, while we reach $N=12$ in the BBFG
method using a mid-to-high-end laptop (DELL XPS 15 with an Intel Core
i7-8750H). If matrix size of $N=20$ is the target, the BBFG method can perform
a full $P$-matrix determination in $\sim 711$ seconds, while it would take
$\sim 15.3\times 10^{6}$ years in the factorial approach. The BBFG method is
the fastest completely general solution for the problem of computing a
$P$-matrix from any $W$-matrix. For several algorithms, the $W$-matrix has
special characteristics that can be exploited to further increase efficiency.
For instance, if by shuffling the rows and columns the $W$-matrix can be made
into a block form, where squared blocks at the diagonal have only zero’s at
their right and upper side, the permanent is equal to the product of the
block’s permanents. For instance, if $W_{14}=W_{24}=W_{34}=0$ we have two
blocks, $3\times 3$ and $1\times 1$. If $W_{13}=W_{14}=W_{23}=W_{24}=0$, we
can identify 2 blocks of $2\times 2$ etc. Identification of blocks can hugely
decrease the computation of a large permanent. Another speed-up can be made if
all rows in the $W$-matrix are a sequence of ones followed by all zeros, or
can be made into that form after previously mentioned column and row
operations. This makes an order ${\mathcal{O}}(n^{2})$ approach possible. We
will further discuss this in Sec. Application: $\infty$RETIS.

The infinite swapping approach changes the aforementioned algorithm from step
3:

  1. 3.

If the worker is finished with its MC move in a specific ensemble, the new
state is accepted or rejected (but not yet sampled) according to Eq. 1. The
ensemble is free.

  2. 4.

Determine the $W$-matrix based on all unoccupied ensembles, calculate the
$P$-matrix based on Eq. 5, and update all the unoccupied ensembles by sampling
all free states with the fractional probabilities corresponding to the columns
in the $P$-matrix.

  3. 5.

Pick randomly one of the free ensembles $e_{j}$.

  4. 6.

Pick one of the available states $(s_{1},s_{2},\ldots)$ based on a weighted
random selection in which state $s_{i}$ has a probability of $P_{ij}$ to be
selected.

  5. 7.

The worker is assigned to do a new standard move in ensemble $e_{j}$ based on
previous state $s_{i}$. Go to step 3.

## Application: $\infty$RETIS

Replica Exchange Transition interface sampling (RETIS) [19, 20] is a
quantitative path sampling algorithm in which the sampled states are short
molecular trajectories (paths) with certain start- and end-conditions, and a
minimal progress condition. New paths are being generated by a Monte Carlo
move in path space, such as the shooting move [21] in which a randomly
selected phase point of the previous path is randomly modified and then
integrated backward and forward in time by means of molecular dynamics (MD).
The required minimal progress increases with the rank of the ensemble such
that the final ensemble contains a reasonable fraction of transition
trajectories. The start- and end-conditions, as well as the minimal progress,
are administered by the crossings of interfaces
$(\lambda_{0},\lambda_{1},\ldots,\lambda_{M})$ with
$\lambda_{k+1}>\lambda_{k}$, that can be viewed as non-intersecting
hypersurfaces in phase space having a fixed value of the reaction coordinate.
A MC move that generates a trial path not fulfilling the path ensemble’s
criteria is always rejected. RETIS defines different path ensembles based on
the the direction of the paths and the interface that has to be crossed, but
all paths start by crossing $\lambda_{0}$ (near the reactant state/state $A$)
and they end by either crossing $\lambda_{0}$ again or reaching the last
interface $\lambda_{M}$ (near the product state/state $B$). There is one
special path ensemble, called $[0^{-}]$, that explores the left side of
$\lambda_{0}$, the reactant well, while all other path ensembles, called
$[k^{+}]$ with $k=0,1,\ldots M-1$, start by moving to the right from
$\lambda_{0}$ reaching at least $\lambda_{k}$.

A central concept in RETIS is the so-called overall crossing probability, the
chance that a path that crosses $\lambda_{0}$ in the positive direction
reaches $\lambda_{M}$ without recrossing $\lambda_{0}$. It provides the rate
of the process when multiplied with the flux through $\lambda_{0}$ (obtained
from the path lengths in $[0^{-}]$ and $[0^{+}]$ [20]) and is usually an
extremely small number. The chance that any of the sampled paths in the
$[0^{+}]$ path ensemble crosses $\lambda_{M}$ is generally negligible, but a
decent fraction of those ($\sim 0.1-0.5$) will cross $\lambda_{1}$ and some
even $\lambda_{2}$. Likewise, paths in the $[k^{+}]$, $k>0$, path ensembles
have a much higher chance to cross $\lambda_{k+1}$ than a $[0^{+}]$-path as
they already cross $\lambda_{k}$. This leads to the calculation of $M$ local
conditional crossing probabilities, the chance to cross $\lambda_{k+1}$ given
$\lambda_{k}$ was crossed for $k=0,1\ldots M-1$, whose product gives an exact
expression for the overall crossing probability with an exponentially reduced
CPU cost compared to MD.

The efficiency is further hugely improved by executing replica exchange moves
between the path ensembles. These swaps are essentially cost-free since there
is no need to simulate additional ensembles that are not already required. An
accepted swapping move in RETIS provides new paths in two ensembles without
the expense of having to do MD steps. The enhancement in efficiency is
generally even larger than one would expect based on those arguments alone as
path ensembles higher-up the barrier provide a similar effect as the high-
temperature ensembles in parallel tempering. In addition, point-exchange moves
between the $[0^{-}]$ and $[0^{+}]$ are performed by exchanging the end- and
start-points of these path that are then continued by MD at the opposite site
of the $\lambda_{0}$ interface.

While TIS [22] (without replica exchange) can run all path ensembles
embarrassingly parallel, the RETIS algorithm increases the CPU-time
efficiency, but is difficult to parallelize and open source path-sampling
codes, like OpenPathSampling [23] and PyRETIS [24], implement RETIS as a fully
sequential algorithm. The path length distributions are generally broad with
an increasing average path length as function of the ensemble’s rank. This
becomes increasingly problematic the more ensembles you have as they all have
to wait for the slowest ensemble. This means that while RETIS will give you
the best statistics per CPU-hour, it might not give you the best statistics in
wall-time. With the continuous increase in computing power, trading some CPU-
time efficiency for wall-time efficiency, getting the answer faster while
spending more CPU-cycles, might be preferential. Our parallel scheme can
effectively deal with the unequal CPU cost of the replicas, which allows us to
increase the wall-time efficiency with no or minimal reduction in CPU-time
efficiency.

The $W$-matrix in RETIS. If there are $M+1$ interfaces,
$\lambda_{0},\lambda_{1},\ldots,\lambda_{M}$, there are also $N=M+1$
ensembles, $[0^{-}],[0^{+}],[1^{+}],\ldots,[(M-1)^{+}]$. For $K$ workers, the
size of the $W$-matrix is, hence, either $(N-K+1)\times(N-K+1)$ or
$(N-K)\times(N-K)$ as swappings are executed when 1 of the $K$ workers is
free, while the remaining $K-1$ workers occupy path ensembles that are locked
and do not participate in the swap. The smallest matrix occurs when one worker
is occupying both $[0^{-}]$ and $[0^{+}]$ during the point exchange move, as
described in the simulation methods.

Paths can be represented by a sequence of time slices, the phase points
visited by the MD trajectory. For a path of length $L+1$,
$X=(x_{0},x_{1},\ldots,x_{L})$, the plain path probability density $\rho(X)$
is given by the probability of the initial phase point times the dynamical
transition probabilities to go from one phase point to the next:
$\rho(X)=\rho(x_{0})\phi(x_{0}\rightarrow x_{1})\phi(x_{1}\rightarrow
x_{2})\ldots\phi(x_{L-1}\rightarrow x_{L})$. Here, the transition
probabilities depend on the type of dynamics (deterministic, Langevin, Nosé-
Hoover dynamics, etc). The weight of a path within a specific path ensemble
$\rho_{j}(X)$ can be expressed as the plain path density times the indicator
function ${\bf 1}_{e_{j}}$ and possibly an additional weight function
$w_{j}(X)$: $\rho_{j}(X)=\rho(X)\times{\bf 1}_{e_{j}}(X)\times w_{j}(X)$. The
indicator function equals 1 if the path $X$ belongs to ensemble $e_{j}$.
Otherwise it is 0. The additional weight function $w_{j}(X)$ is part of the
high-acceptance protocol that is used in combination with the more recent path
generation MC moves such as stone skipping [25] and wire-fencing [26]. Using
these ”high-acceptance weights”, nearly all the CPU intensive moves can be
accepted as they are tuned to cancel the $P_{\rm gen}$-terms in Metropolis-
Hastings scheme, Eq. 1, and the effect of the non-physical weights is undone
in the analysis by weighting each sampled path with the inverse of $w_{j}(X)$.

While the path probability $\rho(s_{i}=X)$ is difficult to compute,
determining ${\bf 1}_{j}(s_{i})$ and $w_{j}(s_{i})$ is trivial. It is
therefore a fortunate coincidence that we can replace $W_{ij}=\rho_{j}(s_{i})$
with

$\displaystyle W_{ij}={\bf 1}_{e_{j}}(s_{i})w_{j}(s_{i})$ (6)

because the $P$-matrix does not change if we divide or multiply a row by the
same number, as mentioned in Sec. Methods. Except for $[0^{-}]$, all path
ensembles have the same start and end condition and only differ with respect
to the interface crossing condition. A path that crosses interface
$\lambda_{k}$ automatically crosses all lower interfaces $\lambda_{l<k}$.
Reversely, if the path does not cross $\lambda_{k}$, it won’t cross any of the
higher interfaces $\lambda_{l>k}$. This implies that if the columns of
$W_{ij}$ are ordered such that the 1st column ($e_{1}$) is the first available
ensemble from the sequence $([0^{-}],[0^{+}],[1^{+}],\ldots,[(M-1)^{+}])$, the
2nd column ($e_{2}$) is the second available ensemble etc, most rows will end
with a series of zeros.

Reordering the rows with respect to the number of trailing zeros, almost
always ensures that the $W$-matrix can be brought into a block-form such that
the permanent can be computed faster based on smaller matrices. In particular,
if $[0^{-}]$ is part of the free ensembles, it will always form a $1\times 1$
block as there is always one and no more than one available path that fits in
this ensemble.

If high-acceptance is not applied, we have $w_{j}(X)=1$ and each row in the
$W$-matrix (after separating the $[0^{-}]$ ensemble if it is part of the free
ensembles) is a sequence of ones followed by all zeros. The $W$-matrix can
hence be represented by an array $(n_{1},n_{2},n_{3},\ldots n_{n})$ where each
integer $n_{i}$ indicates the number of ones in row $i$. As we show in the SI,
the permanent of such a $W$-matrix is simply the product of $(n_{i}+1-i)$:
${\rm perm}(W)=\prod_{i}(n_{i}+1-i)$. Further, the $P$-matrix can be
constructed from following order ${\mathcal{O}}(n^{2})$ method.

The first step is to order the rows of the $W$-matrix such that $n_{1}\leq
n_{2}\leq\ldots\leq n_{n}$. We then fill in the $P$-matrix from top to bottom
for each row using

$\displaystyle P_{ij}=\begin{cases}0,\text{ if }W_{ij}=0\\\
\frac{1}{n_{i}+1-i},\text{ if }W_{ij}=1\text{ and }[W_{(i-1)j}=0\text{ or
}i=1]\\\ \left(\frac{n_{i-1}+1-i}{n_{i}-i}\right)P_{(i-1)j},\text{ otherwise
}\end{cases}$ (7)

The approach is extremely fast and allows the computation of $P$-matrices from
a large $W$-matrix, up to several thousands, within a second of CPU-time. The
above method applies whenever the rows of the $W$-matrix can be transformed
into sequence of ones followed by all zeros. Besides RETIS without high-
acceptance, this would apply to other MC methods like subset-sampling [27] or
umbrella sampling [28] with semi-infinite rectangular windows.

## Results and Discussion

To test our algorithms we ran 3 types of simulations. First a memoryless
single variable stochastic (MSVS) process was simulated in order to mimic a
RETIS simulation in which the average path length increases linearly with the
rank of the ensemble. A ”path” is created by drawing 2 random numbers where
the first determines how much progress a path makes and the second determines
the path length. These two outcomes are variable and depend on the rank of the
ensemble such that the fictitious path in ensemble $[k^{+}]$ has a 0.1
probability to cross $\lambda_{k+1}$ and has an average path length of
approximately $k/10$ seconds (see Section Materials and Methods). The worker
is paused for a number of seconds equal to the path length before it can
participate in replica moves to mimic the time it would take to do all the
necessary MD steps. While this artificial simulation allows us to investigate
the potential strength of the method to tackle extremely rare events, it
cannot reveal the effect of correlations between accepted paths when fast
exploration of the reaction coordinate’s orthogonal directions are crucial. To
analyze this effect, we also ran a 2D membrane permeation system with two
slightly asymmetric channels [29]. Lastly, to study our algorithm with a more
generic $W$-matrix that needs to be solved via BBFG formula, we also ran a set
of underdamped Langevin simulations of a particle in a double well potential
[30] using the recent wire fencing algorithm with the high acceptance protocol
[26]. All simulation results were performed using 5 independent runs of 12
hours. Errors were based on the standard deviations from these 5 simulations,
except for the MSVS process, where a more reliable statistical error was
desired for the comparison with analytical results. Here, block errors were
determined on each of the five simulations based on the running average of the
overall crossing probability. The block errors were finally combined to obtain
the statistical error in the average of the five simulations.

### Memoryless single variable stochastic (MSVS) process

Table 1 reports the overall crossing probabilities and their statistical
errors for a system with 50 interfaces and 1, 5, 10, 15, 20, 25, 30, 35, 40,
45, and 50 workers.

Table 1: Results of the 3 model systems showing crossing probabilities ($P_{\rm cross}$), permeabilities (perm.), and rates for different number of workers (#w). All results are shown in dimensionless units. Errors are based on single standard deviations. Values shown in the lower part are a: exact result, b: Ref. [29], c: approximated value based on Kramers’ therory (see SI), d: Ref. [30], and e: Ref. [26].  MSVS | two-channel system | double well with wire fencing
---|---|---
#w | $P_{\rm cross}/10^{-50}$ | #w | $P_{\rm cross}/10^{-5}$ | perm./$10^{-6}$ | #w | $P_{\rm cross}/10^{-7}$ | rate/$10^{-7}$
1 | $0.61\pm 0.33$ | 1 | $1.52\pm 0.17$ | $1.28\pm 0.14$ | 1 | $5.91\pm 0.18$ | $2.59\pm 0.07$
5 | $1.47\pm 1.04$ | 2 | $1.63\pm 0.24$ | $1.37\pm 0.20$ | 2 | $5.70\pm 0.13$ | $2.51\pm 0.06$
10 | $0.86\pm 0.51$ | 3 | $1.52\pm 0.07$ | $1.28\pm 0.06$ | 3 | $5.57\pm 0.19$ | $2.45\pm 0.08$
15 | $0.68\pm 0.08$ | 4 | $1.42\pm 0.10$ | $1.19\pm 0.08$ | 4 | $5.20\pm 0.30$ | $2.34\pm 0.12$
20 | $1.02\pm 0.13$ | 5 | $1.40\pm 0.12$ | $1.18\pm 0.10$ | 5 | $5.05\pm 0.41$ | $2.23\pm 0.18$
25 | $1.02\pm 0.17$ | 6 | $1.54\pm 0.06$ | $1.30\pm 0.05$ | 6 | $5.49\pm 0.29$ | $2.42\pm 0.13$
30 | $1.26\pm 0.24$ | 7 | $1.48\pm 0.08$ | $1.24\pm 0.07$ | 7 | $4.99\pm 0.39$ | $2.21\pm 0.17$
35 | $1.05\pm 0.15$ | 8 | $1.46\pm 0.08$ | $1.23\pm 0.06$ | 8 | $4.88\pm 0.43$ | $2.15\pm 0.19$
40 | $1.05\pm 0.14$ | 9 | $1.42\pm 0.10$ | $1.20\pm 0.08$ |  |  |
45 | $0.93\pm 0.09$ | 10 | $1.44\pm 0.08$ | $1.21\pm 0.07$ |  |  |
50 | $1.00\pm 0.07$ | 11 | $1.41\pm 0.09$ | $1.19\pm 0.08$ |  |  |
|  | 12 | $1.30\pm 0.15$ | $1.09\pm 0.12$ |  |  |
literature/theoretical result
|  |  | $1.23\pm 0.16^{\rm b}$ | $1.06\pm 0.14^{\rm b}$ |  |  | $2.79\pm 0.70^{\rm d}$
|  |  |  |  |  | $5.84\pm 0.13^{\rm e}$ | $2.58\pm 0.06^{\rm e}$
| $1.00^{\rm a}$ |  | $1.61^{\rm c}$ | $1.37^{\rm c}$ |  | $5.83^{\rm c}$ | $2.58^{\rm c}$

All values are within a 50% deviation from the true value of $10^{-50}$ with
the more accurate estimates for the simulations having a large number of
workers. Also, the true value is within one standard deviation of the reported
averages for 70% of the data points, as is expected from the standard Gaussian
confidence intervals. Figure 1a shows the scaling of the MD time (solid lines)
and number of MC moves (dashed lines) of the MSVS simulations (orange)
compared to linear scaling (black) and the expected scaling for standard
replica exchange (REPEX) in which ensembles are updated in cohort (purple).

Whereas the number of ”MD steps” and MC moves quickly levels off to a nearly
flat plateau in the standard approach due to workers being idle as they need
to wait for the slowest worker, the replica exchange approach developed in
this article shows a perfect linear scaling with respect to the MD time. The
number of MC moves in the new method shows an even better than linear scaling
due to the fact that the ensembles with shorter ”path lengths” get simulated
relatively more often with more workers, resulting in more MC moves per
second. This in itself does not necessarily mean that the simulations converge
much faster because the additional computational effort may not be targeted to
the sampling where it is needed. If we neglect the fact that path ensemble
simulations are correlated via the replica exchange moves, we can write that
the relative error in overall crossing probability $\epsilon$ follows from the
relative errors in each path ensemble $\epsilon_{i}$ via:
$\epsilon^{2}=\sum_{i}\epsilon_{i}^{2}$. It is henceforth clear that
additional computational power should not aim to lower the error in a few path
ensembles that were already low compared to other path ensembles. We therefore
measure the effectiveness of the additional workers by calculating
computational efficiencies. The efficiency of a specific computational method
is here defined as the inverse computer time, CPU- or wall-time, to obtain an
overall relative error equal to 1: $\epsilon=1$.

In figure 1d the efficiencies based on wall-time (solid) and CPU-time (dashed)
are plotted for the MSVS process. These plots depends on the ability of
computing reliable statistical errors in the overall crossing probability that
is an extremely small number, $10^{-50}$. The somewhat fluctuating behavior of
these curves should hence be viewed as statistical noise as the confidence
interval of these efficiencies depends on the statistical error of this error.
Despite that, clear trends can be observed in which the CPU-time efficiency is
more or less flat, while the wall-time efficiency shows an upward trend. If we
neglect the effect of replica exchange moves on the efficiency, we can relate
these numerical results with theoretical ones [31, 20] for any possible
division of a fixed total CPU-time over the different ensembles. A common
sense approach would be to aim for the same error $\epsilon_{i}$ in each
ensemble (which implies doing the same number of MC moves per ensemble) or to
divide the total CPU-time evenly over the ensembles. These two strategies
correspond to the case $K=1$ or standard RETIS and $K=N$ or standard TIS,
respectively. Ref. [31] showed that these two strategies provide the same
efficiency and in the SI we derive that this leads to a wall-time efficiency
as function of the number of workers ($K$) equal to $K/56250$ which is the
continuous purple line in figure 1d. The optimum division, however, would give
a slightly better wall-time efficiency equal to $K/50000$ which is the
continuous black line in this figure. Also shown in figure 1d are the expected
theoretical efficiencies based on the numerical distribution of MC moves in
each ensemble. This hybrid numerical/theoretical result is shown by the small
purple dots. This shows that $\infty$RETIS, at least for a system in which the
path length grows linearly with the ensemble’s rank, naturally provides a
division of the computational resources that is even better than TIS ($K=N$)
or RETIS $(K=1)$. Yet, due to statistical inaccuracies this is only evident
for the $K=15$ case. The best wall-time efficiency is obtained for the case
$K=N$, which is essentially equivalent of running independent TIS simulations
(i.e. without doing any replica exchange moves). We do not expect this to
apply to more complex systems where the replica exchange move is a proven
weapon for efficient sampling.

### Two-channel simulations

In the middle column of table 1 we report the calculated crossing
probabilities and permeabilities for 5 simulations for every number of
workers. All simulations are somewhat higher, though still in good agreement
with the previous simulation from Ref. [29]. We also evaluated the approximate
result based on Kramers’ theory (see SI) which seem to confirm the results
obtained in this paper.

Figure 1b shows the scaling of the MD time (solid lines) and number of MC
moves (dashed lines) of the two-channel simulations (blue) compared to linear
scaling (black). We see a slightly worse than linear scaling of the MD time,
which might just be due to a small positive fluctuation of the 1 worker data-
point. We also see a similar more than linear scaling in the number of MC
moves as with the MSVS simulations, for the same reason. In figure 1e the
efficiencies based on wall-time (solid) and CPU-time (dashed) are plotted for
the two-channel system. The CPU-time efficiency is more or less flat until 8
workers after which it starts to drop off. The wall-time efficiency shows an
upward trend until 10 workers after which it starts to drop off as well. We
assign this drop to the reduction of replica exchange moves which is an
essential aspect for sampling this system efficiently [29]. This is tangible
from figure S1 in the SI where we plot fraction of trajectories, passing
through $\lambda_{M-1}$, that are in the lower barrier channel. While from the
average fraction it still looks like the simulations sampled both channels for
any number of workers, 4 out of the 5 simulations in the $K=N=12$ case solely
visited one of the two channels. This is in agreement with previous TIS
results [29]. The $K=11$ case already provides a dramatic improvement, but is
still expected to be sub optimal due to the relatively low frequency of
replica exchange moves compared to $K<11$. As reported in ref. [29], this
ratio requires many MC moves to converge to the theoretical value of 0.71
without the added MC moves introduced in that paper. We did not simulate with
these added moves and thus see the same slow convergence for all of our
simulations. From this 2D system it would indicate that having half the number
of ensembles as workers is a safe bet for optimum efficiency.

a) b) c)
d) e) f)

Figure 1: The average scaling of total MD time (cumulative time spend by all
the workers) (solid) and MC moves (dashed) (a-c) and wall-time (solid) and the
CPU-time (dashed) efficiencies (d-f) for each number of workers. This is shown
for the memoryless single variable stochastic (MSVS) process (a, d, orange),
the two-channel system (b, e, blue) and the double well with wire fencing (c,
f, green) simulations. Each of the data points is based on 5 independent
simulations. For the scaling plots, the black lines are guides for linear
scaling from the 1 worker data-point. The purple lines in the scaling plot for
the MSVS simulations (a) show what the scaling would be if we had to wait for
the slowest ensemble to finish for each MC move. The black line, purple line,
and points in the efficiency plot of the MSVS process (d) show the optimal,
optimal TIS/RETIS, and hybrid wall-time efficiency, respectively, as computed
in the SI.

### Double well 1D barrier using wire fencing

In the right column of table 1 we report the calculated crossing probabilities
and rates for the underdamped Langevin particle in the 1D double well
potential. All simulations are in reasonable agreement with each other and the
results of Refs. [26] and and [30], as well as the approximate value based on
Kramers’ theory. However, while these results confirm the soundness of the
method, the scaling and efficiency are less convincing. Figure 1b shows a
significantly worse than linear scaling. On further inspection we found the
average time per MC move was significantly smaller than our infinite-swapping
goal of 1 s when the simulation was run with more than 2 workers. This results
in a bottleneck on how many MC moves can be started per second, which is the
reason for the observed bad scaling. It still is slightly positive instead of
flat as the infinite swapping procedure becomes quicker with more workers due
to the smaller $W$-matrix. The same bottleneck can be seen in figure 1f were
both efficiencies plummet with more than 2 workers. The reported scaling
deficiency is of little significance for actual molecular systems where the
creation of a full path takes minutes to hours rather than subseconds.

## Conclusions

We developed a new generic replica exchange method that is able to effectively
deal with MC moves with varying CPU costs, for instance due to the algorithmic
complexity of the MC moves. An essential aspect of the method is that the
number of workers, who execute the ensemble’s specific MC moves in parallel,
is less than the number of ensembles. Once a worker is finished with its move,
replica exchange moves are carried out solely between those ensembles that are
not occupied by a worker. This implies that the ensembles are updated at
irregular intervals and a different number of MC moves will be executed for
each ensemble. As a result, the conceptual viewpoint in which the set of
replica’s are viewed as a single superstate is no longer valid and the
existence of some kind of detailed-balance relation is no longer trivial. To
prove the exactness of our approach, we introduced some new conceptual views
on the replica exchange methodology that is different from the common
superstate principle. Instead, we show that the distributions in the new
approach are conserved for each ensemble individually via a twisted detailed
balance relation in which the other ensembles constitute an environment that
is potentially actively involved in the MC move of the ensemble considered. In
addition, the method can be combined with an infinite swapping approach
without the factorial scaling based on a mathematical reformulation using
permanents.

We applied the novel replica exchange technique on a path sampling algorithm,
RETIS, which is a prototype of algorithm where the costs for a Monte Carlo
move can vary enormously. The resulting new path sampling algorithm, coined
$\infty$RETIS, was thereafter tested on three model systems. The results of
these simulations show that the number of MD steps increase linearly with the
number of workers invoked as long as the ensemble’s MC move has a lower
computational cost than the replica exchange move carried out by the
scheduler. The number of executed MC moves shows an even better than linear
scaling. Moreover, the efficiency increases linearly with the number of
workers for a low-dimensional system in which the replica exchange has little
effect, while it has an optimum in more complex systems as the number of
successful replica exchange moves decreases when the number of workers is
close to the number of ensembles.

In summary, the replica exchange method discussed in this paper has a clear
potential to accelerate present path sampling simulations, but can also be
combined with many other complex algorithms including those that are yet to be
invented. With continuing trend to to run progressively more massively-
parallel computing jobs, our algorithm is likely to gain importance and will
open up new avenues in the field of molecular simulations and beyond.

### Supporting Information Appendix (SI)

This article contains Supplementary Information

## Materials and Methods

### Simulation methods

The implementation of $\infty$RETIS was structured as follows. We start $1\leq
K\leq N$ worker- and 1 scheduler-process. Each of the worker-processes is
going to process ensemble specific MC moves while the scheduler-process will
do all the replica exchange moves and submits new jobs to the workers. All
ensembles/trajectories that are currently being updated by a worker-process
are not considered for MC moves by the scheduler, essentially being ’locked’.
This means that no data is written for those ensembles and they are not valid
targets for swapping moves. After a worker is done, it submits the result to
the scheduler, the scheduler then unlocks the returned ensemble/trajectory and
executes the replica exchange moves on all ensembles/trajectories that that
are not locked. It then submits a new job to the freed worker for performing a
new MC move in a randomly chosen free ensemble (or two ensembles in case of a
point exchange move) and locks the involved ensembles/trajectories.

In the $\infty$RETIS method there are two kind of ensemble moves that involve
MD steps. The first one is the shooting move (either standard shooting [21] or
the more recent sub-trajectory moves [25, 26]) in which a new path is being
generated from an old path within a single ensemble. The second one is the
point exchange move between $[0^{-}]$ and $[0^{+}]$. If a worker is assigned
to this task, it means that both $[0^{-}]$ and $[0^{+}]$ are occupied by this
worker. The scheduler ensures that there is never more than 1 worker
considered free at a given time. When the free worker is assigned to perform a
new MC move, each of the ensembles have an equal probability to be selected.
If $[0^{+}]$ or $[0^{-}]$ is selected and the other is also free, there is a
50% chance to perform a $[0^{-}]\leftrightarrow[0^{+}]$ point exchange move
instead of a shooting move in the selected ensemble.

### Memoryless single variable stochastic (MSVS) process

No actual MD is run for the MSVS simulations. Instead, we directly sample two
random numbers, $r_{1}$ and $r_{2}$ from an uniform distribution $\in[0,1)$ to
set the path’s progress and the path length. A path in ensemble $[k^{+}]$ is
assumed to cross interface $\lambda_{k+l}$ if $r_{1}<(0.1)^{l}$. After this,
we wait a random time, $t=0.2\,r_{2}\,k+0.1$ in seconds. This was done to
simulate both the increasing average simulation time and variance for outer
ensembles. This setup means that we have no history dependence and allows us
to compute the theoretical values show in figure 1. 5 independent
$\infty$RETIS simulations were run with $1,5,10,15,\ldots,45,50$ workers.

### Double channel simulations

In order to investigate the effect of our algorithm on the ergodicity of the
sampling, a 2D two-channel simulation was run as described in reference [29].
The new RETIS moves introduced in that paper (mirror-move and target-swap
move) were not used. Instead, MD was only run to do shooting moves or the
$[0^{-}]\leftrightarrow[0^{+}]$ point exchanges. As the MD for this system
completed too fast, every worker was set to wait $9$ times the time it took to
run the MD before returning the result. 5 independent $\infty$RETIS
simulations were run with $1,2,\ldots,11,12$ workers.

### 1D double well with wire fencing

In order to investigate the accuracy with a $W$ matrix that contains more
numbers than $0$s or $1$s we simulated a 1D double-well system [30] together
with the high-acceptance version of a novel path-sampling algorithm, wire
fencing. The algorithm is described in reference [26], but for us the relevant
part is that the high-acceptance weight is the number of frames that a path
has outside the interface for each ensemble times an extra factor 2 if the
path ends at the last interface. As for the two-channel system, a worker was
set to wait 9 times the time it took to complete the MD move before returning
the result. 5 independent $\infty$RETIS simulations were run with
$1,2,\ldots,7,8$ workers with interfaces placed at
$[-0.99,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,1.0]$.

## ACKNOWLEDGMENTS

We acknowledge funding from the Research Council of Norway through FRINATEK
Project No. 275506.

#### The authors declare no conflict of interest

## References

  * Metropolis  _et al._ [1953] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21, 1087 (1953).
  * Hastings [1970] W. Hastings, Monte-Carlo sampling methods using Markov chains and their applications, Biometrika 57, 97 (1970).
  * Swendsen and Wang [1986] R. H. Swendsen and J. S. Wang, Replica monte-carlo simulation of spin-glasses, Phys. Rev. Lett. 57, 2607 (1986).
  * Marinari and Parisi [1992] E. Marinari and G. Parisi, Simulated tempering - a new monte-carlo scheme, Europhysics Lett. 19, 451 (1992).
  * Sugita and Okamoto [1999] Y. Sugita and Y. Okamoto, Replica-exchange molecular dynamics method for protein folding, Chem. Phys. Lett. 314, 141 (1999).
  * Siepmann and Frenkel [1992] J. I. Siepmann and D. Frenkel, Configurational bias Monte-Carlo - a new sampling scheme for flexible chains, Molecular Physics 75, 59 (1992).
  * Vlugt  _et al._ [1999] T. Vlugt, R. Krishna, and B. Smit, Molecular simulations of adsorption isotherms for linear and branched alkanes and their mixtures in silicalite, J. Phys. Chem. B 103, 1102 (1999).
  * Frenkel and Smit [2002] D. Frenkel and B. Smit,  _Understanding molecular simulations from algorithms to applications_ (Academic press, San Diego, California, U.S.A., 2002).
  * Dellago  _et al._ [1998a] C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, Transition path sampling and the calculation of rate constants, J. Chem. Phys. 108, 1964 (1998a).
  * Swendsen and Wang [1987] R. H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in monte carlo simulations, Phys. Rev. Lett. 58, 86 (1987).
  * Peters and de With [2012] E. A. J. F. Peters and G. de With, Rejection-free monte carlo sampling for general potentials, Phys. Rev. E 85, 026703 (2012).
  * Michel  _et al._ [2014] M. Michel, S. C. Kapfer, and W. Krauth, Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps, J. Chem. Phys. 140, 054116 (2014).
  * Plattner  _et al._ [2011] N. Plattner, J. D. Doll, P. Dupuis, H. Wang, Y. Liu, and J. E. Gubernatis, An infinite swapping approach to the rare-event sampling problem, The Journal of Chemical Physics 135, 134111 (2011), https://doi.org/10.1063/1.3643325 .
  * Plattner  _et al._ [2013] N. Plattner, J. D. Doll, and M. Meuwly, Overcoming the rare event sampling problem in biological systems with infinite swapping, J. Chem. Theory Comput. 9, 4215 (2013).
  * Lu and Vanden-Eijnden [2019] J. Lu and E. Vanden-Eijnden, Methodological and computational aspects of parallel tempering methods in the infinite swapping limit, J Stat Phys 174, 715 (2019).
  * Balasubramanian [1980] K. Balasubramanian,  _Combinatorics and diagonals of matrices_ , Ph.D. thesis, Loyola College, Madras, India (1980).
  * Bax [1998] E. Bax,  _Finite-difference Algorithms for Counting Problems_ , Ph.D. thesis, California Institute of Technology, Pasadena, United States of America (1998).
  * Glynn [2010] D. G. Glynn, The permanent of a square matrix, European Journal of Combinatorics 31, 1887 (2010).
  * van Erp [2007] T. van Erp, Reaction rate calculation by parallel path swapping, Phys. Rev. Lett. 98, 268301 (2007).
  * Cabriolu  _et al._ [2017] R. Cabriolu, K. M. S. Refsnes, P. G. Bolhuis, and T. S. van Erp, Foundations and latest advances in replica exchange transition interface sampling, J. Chem. Phys. 147, 152722 (2017).
  * Dellago  _et al._ [1998b] C. Dellago, P. G. Bolhuis, and D. Chandler, Efficient transition path sampling: Application to lennard-jones cluster rearrangements, The Journal of Chemical Physics 108, 9236 (1998b), https://doi.org/10.1063/1.476378 .
  * van Erp  _et al._ [2003] T. S. van Erp, D. Moroni, and P. G. Bolhuis, A novel path sampling method for the sampling of rate constants, J. Chem. Phys. 118, 7762 (2003).
  * Swenson  _et al._ [2019] D. W. H. Swenson, J.-H. Prinz, F. Noe, J. D. Chodera, and P. G. Bolhuis, Openpathsampling: A python framework for path sampling simulations. 2. building and customizing path ensembles and sample schemes, J. Chem. Theory Comput. 15, 837 (2019).
  * Riccardi  _et al._ [2020] E. Riccardi, A. Lervik, S. Roet, O. Aaroen, and T. S. van Erp, Pyretis 2: An improbability drive for rare events, J. Comput. Chem. 41, 370 (2020).
  * Riccardi  _et al._ [2017] E. Riccardi, O. Dahlen, and T. S. van Erp, Fast decorrelating monte carlo moves for efficient path sampling, J. Phys. Chem. Lett. 8, 4456 (2017).
  * Zhang  _et al._ [2021] D. T. Zhang, E. Riccardi, and T. S. van Erp, Path sampling with sub-trajectory moves, In preparation xx, xx (2021).
  * Au and Beck [2001] S.-K. Au and J. L. Beck, Estimation of small failure probabilities in high dimensions by subset simulation, Probabilistic Eng. Mech. 16, 263 (2001).
  * Torrie and Valleau [1977] G. Torrie and J. Valleau, Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella samping, J. Comp. Phys. 23, 187 (1977).
  * Ghysels  _et al._ [2021] A. Ghysels, S. Roet, S. Davoudi, and T. S. van Erp, Exact non-markovian permeability from rare event simulations, Phys. Rev. Research 3, 033068 (2021).
  * van Erp [2012] T. van Erp, Dynamical rare event simulation techniques for equilibrium and nonequilibrium systems, Adv. Chem. Phys. 151, 27 (2012).
  * van Erp [2006] T. S. van Erp, Efficiency analysis of reaction rate calculation methods using analytical models I: The two-dimensional sharp barrier, J. Chem. Phys. 125, 174106 (2006).