extfieldsofstudy: Medicine Physics

text:
Supervised and Unsupervised Machine Learning of Structural Phases of Polymers Adsorbed to Nanowires

We identify configurational phases and structural transitions in a polymer nanotube composite by means of machine learning. We employ various unsupervised dimensionality reduction methods, conventional neural networks, as well as the confusion method, an unsupervised neural-network-based approach. We find neural networks are able to reliably recognize all configurational phases that have been found previously in experiment and simulation. Furthermore, we locate the boundaries between configurational phases in a way that removes human intuition or bias. This could be done before only by relying on preconceived, ad-hoc order parameters.

I. INTRODUCTION
We previously studied soft-solid matter nano-composites, in particular systems composed of flexible polymers adsorbed at thin nanostrings or tubes [1][2][3]. Such systems are believed to play an important role in current and future development of high-performance nanomaterials. Carbon nanotubes, for example, have been functionalized by wrapping them with certain types of polymers to serve as biosensors for the detection of glucose [4,5]. However, the successful fabrication of such materials depends on a variety of parameters. In particular, the wetting behavior of carbon nanotubes has been shown to be one critical parameter for the development of nanocomposites [6]. We have previously developed and employed a coarse-grained model to investigate nanoscale wetting and adhesion phenomena using Monte Carlo methods; we identified various, structurally different low-temperature phases including globular polymers simply attached to the nanostring and polymers completely wrapping, or coating, the substrate [1]. One particular problem that we recognized was the classification of structural phases at low temperatures depending on various model parameters. We have addressed the problem in the past by the ad hoc introduction of order parameters to identify boundaries between such structural phases. In this paper we show how a less biased approach, based on machine learning, can be deployed. We revisit the earlier introduced configurational phase diagram, aiming at identifying classes of the polymer-wire system and the boundaries between them without any assumptions or other input based on * Current address: School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, U.S.A.; Correspondence email address: qparker3@gatech.edu a human perception of structure. In a more general context, such automated structure identification can also provide means to recognize system configurations during Monte Carlo sampling. This could prove beneficial in order to run generalized-ensemble simulations where different structural phases are assigned individual weights [7], or simply to collect statistics for individual phases during a simulation. Recent years have witnessed significant advances in the use of machine learning (ML) methods for phase classification. In this regard, supervised learning approaches [8][9][10][11][12][13][14], for which the prior labeling of configurations is required, as well as unsupervised learning approaches [15][16][17][18][19][20], which work without such prior labeling, have been attempted. It has been demonstrated, for example, that neural networks (NNs) trained with labeled configurations can encode information about the ordered and disordered phases in model systems by learning the relevant order parameters [8]. In particular, approaches based on purposefully mislabeling configurations and evaluating the network performance have been developed to detect phase transitions [21]. Such a method does not require true labels to be known in advance and therefore no prior knowledge about the existence (or lack thereof) of a transition is needed. It has also been demonstrated to work in the presence of multiple transitions [22].
In the context of unsupervised learning approaches, dimensionality reduction techniques, for example, principal component analysis (PCA), multidimensional scaling (MDS), t-distributed stochastic neighbor embedding (t-SNE), autoencoders, etc., have been found useful in distinguishing ordered and disordered phases [15][16][17][18][19]23]. For systems with clear order parameters, such as the Ising model or the XY model, the latent parameters or the dominant principal components have been shown to directly correlate with the respective order parameters [17,19]. Besides structure recognition, NNs can be trained to predict macroscopic physics quantities such as the total energy, and microscopic quantities such as charge density and magnetization locally for each atom [24]. NNs are also used to learn interatomic potentials, for example. They have been trained to generate effective many-body potentials from ab-initio data [25] and were successfully applied to construct precise phase diagrams of water in molecular-dynamics (MD) simulations over a large range of temperatures and pressures [26]. Another ML potential, ANI-Al, was trained to obtain quantum-level accuracy and has been successfully combined with MD simulations to study shock physics in metals [27]. Training ML surrogate models is also becoming a useful technique to bridge different length and time scales in computer simulations, see [28], for example.
Finally, NNs have been applied in the field of polymer model simulations to study transition between coil and globule structures and recognize Mackay-Anti-Mackay structures, for example [29,30]. Transitions between such crystalline structures in the solid phase are notoriously hard to simulate [31] and advanced generalized ensemble methods have been developed to do so in the past [7,32]. These studies emphasize the benefit of knowing the conformational state of a system during the simulation and ML could contribute valuable information if structures can be reliably recognized. In this paper we will provide more evidence that this can indeed be achieved by employing NNs in the supervised recognition of low-energy configurations of polymers absorbed to a substrate (Sec. III A). Furthermore, we will show how unsupervised ML method can be employed if no previous knowledge of structural phases of a model is available beforehand (Sec. III B). Finally, we will determine boundaries between phases in the model parameter space by training NNs in a conventional way, but also by applying the more recently developed confusion method (Sec. IV).

II. MODEL AND OBSERVED STRUCTURAL PHASES
To model the nanotube-polymer composite we use a coarse-grained bead-spring description for the polymer [33] and an attractive interaction between the monomers and the one-dimensional, continuous string that is derived from a Lennard-Jones potential [1][2][3]. The latter contains two parameters, the effective thickness of the string, σ f , and its attraction strength, f : where r ⊥ is the perpendicular distance between a monomer and the string. The interaction between nonbonded monomers is described by a standard Lennard-Jones potential and there is a weak bending stiffness for consecutive monomer-monomer bonds, as often employed in bead-spring polymer models [33][34][35]. For a more detailed discussion of different approaches to model a thin, cylindrical substrate, see [36]. Depending on the parameter set {σ f , f } in Eqn. (1) low-energy structures take qualitatively different shapes that can be grouped into structural classes or phases. In our previous work we distinguished between four such phases and labeled them Ge, Gi, C, and B; see Fig. 1 for visualisations. Structures like the ones we find in the Gi, Ge, and C regions have been found and imaged in experimental studies before, in particular "clam-shell" (C) polymer nanodroplets have been emphasized [6, Fig. 4]. Globular configurations in phase Ge are similar to structures seen during dewetting of polymers on the surface of carbon nanotubes (CNTs) [6, Fig. 2] while Gi configurations for large values of f show similarities with "barrel-type" nanodroplets [ibid ]. Note that pure monolayer barrel structures (B) can be mapped onto different types of CNTs [37,38]. In fact, we found that region B contains sub-phases with different chiralities corresponding to those found in CNTs [2]. 1

III. MACHINE LEARNING OF STRUCTURAL PHASES
In the following we investigate different supervised and unsupervised machine learning (ML) methods for structure recognition. In ML one typically desires large datasets to reliably train a robust model. However, in the research presented here the data is intrinsically hard to generate since we are analysing states that dominate canonical ensembles at very low temperatures. We use Wang-Landau (WL) sampling [39] to produce these low-energy configurations. Even though WL reliably finds these states, we face the challenge to have to collect many, very different and ideally uncorrelated lowenergy configurations for all parameter values {σ f , f } (see Eq. 1). In an extreme approach and as a proof of concept, we here only record one configuration every time the WL walker explores a low-energy valley and then wait for the walker to move to regions in the phase space corresponding to high temperatures before collecting data again at low energies. Admittedly, such a strategy is computationally expensive and even though applied to generate the dataset analysed in this section, it might not be necessary in that extreme way (see a discussion below in Sec. IV).
In all our simulations, a polymer configuration is represented by the three spatial coordinates of 100 monomers. We use these 300 coordinates (either in raw format or preprocessed, see below) as the feature set for the machine learning algorithms. Although it is possible to utilize an engineered feature set based on our physical intuition of the system (for example by including macroscopic physical observables like the radius of gyration, end-to-end distance, energy, etc.), avoiding such engineered features leaves the machine learning algorithm unbiased and free of any preconceived notions.

A. Supervised Learning: Structure Recognition with Neural Networks
The neural network (NN) is set up with an input layer of 300 neurons, two hidden layers of 50 neurons each and one output layer of four neurons, see Fig. 2. The dataset consists of about 3000 configurations or samples for each polymer type. Two thirds of the dataset are allocated for training, while the remaining data is used for testing. The rectified linear unit (ReLU) activation function is used for the hidden layers and the softmax function for the output layer. We use the Nadam optimization algorithm over ninety epochs for training. Finally, to mitigate overfitting we employ L2 kernel regularization. The results of the NN classification are shown in Fig. 3 where we plot the confusion matrices and the learning curves.
We start by running an analysis with the adsorbed polymer configuration types B, C, Ge, Gi (see Fig. 1) as part of the dataset. We preprocess the data by spatially shifting the monomers in the z direction such that the center of mass of each polymer lies on the xy plane. Figure 3 (a) shows the training and validation accuracy measured at each epoch during the NN training. Both the training accuracy and the validation accuracy rapidly converge to a steady value within 20 epochs. In addition, the training and validation curves are quite close to each other and therefore show no noticeable sign of overfitting. Figure 3 (b) shows the confusion matrix obtained for the validation set, normalized by the number of elements in each class. We see that the offdiagonal elements are zero, except for a very few misclassifications of C-type polymers, yielding an almost 100% overall validation accuracy.
Since the neural network was able to reliably identify all adsorbed polymer structures we also included hightemperature, random-coil polymer structures ("R") not adsorbed to the string as another type and added an output neuron accordingly. Figures 3 (c) and (d) show the corresponding accuracy curves and the confusion matrix. The training and validation accuracies again converge to 1.0, although the convergence rate is slower compared to the previous case. Again, the curves do not show evidence of noticeable overfitting. The slightly increased presence of off-diagonal entries in the "R" row of the confusion matrix indicates a somewhat higher tendency for random coil configurations to be misclassified as other polymer types. However, this is expected as those polymers are random configurations that could, in fact, loosely resemble any of the other classes by chance.

B. Unsupervised Learning: Dimensionality Reduction Methods
Dimensionality reduction methods refer to a class of unsupervised machine-learning techniques that map data from an original, high-dimensional space to a lower-dimensional space while ideally preserving some of the salient properties of the data. In the context of thermodynamic phase classification, for example, such low-dimensional representations of the configuration space have been used to facilitate the visual identification of distinct phases [15,16,19] and to provide insight into the relationship between important features and order parameters of complex systems [17,19,30].
Principal component analysis (PCA) [40] is a linear dimensionality reduction technique which identifies a set of mutually orthogonal unit vectors in a given feature space. These vectors are ordered according to the variance of the data in the corresponding directions, such that the first unit vector indicates the direction of greatest variance in the data. This direction is then called the first principal component, the one with the second highest variance the second principal component, and so on. The principal components are the eigenvectors of the covariance matrix of the data, and hence can be determined by the eigendecomposition of that matrix or the singular value decomposition of the data matrix. The original configurations are then projected into a space spanned by the first m principal components to obtain the desired lower m-dimensional representations. For some spin systems, the principal components have been shown to recover the physical order parameters for phase transitions [17,19].
In addition to PCA, we apply a number of other nonlinear dimensionality reduction methods, namely, multidimensional scaling (MDS) [41], t-distributed stochastic neighbor embedding (t-SNE) [42], Isomap [43], and diffusion map [44]. Note that Isomap becomes equivalent to PCA as the neighborhood size approaches the sample size. Therefore, we limited the neighborhood size to 20 for this demonstration, but also confirmed that changing this number will not change the qualitative results. In general, such non-linear methods identify lower-dimensional manifolds embedded within the higher-dimensional feature space, in which similar data points are clustered together. Typically, manifold learning methods can capture nonlinear relationships within the data that cannot be captured through principal component analysis.

Data Pre-processing
When employing unsupervised learning methods the data typically has to be prepared in some way to obtain meaningful results. In the raw data the polymer is adsorbed at the string at an arbitrary position, while the string is always located at the z-axis in Cartesian coordinates. We here utilize different scaling and coordinate transformation methods to potentially make the features of the polymers more comparable for the machine. A common method in machine learning, referred to as "standard scaling" [45] aims at bringing all features (in our case, monomer coordinates) onto the same length scale by subtracting the mean of all data from each feature and individually scaling each feature to unit variance. Two other ways that do not alter the overall shape of the polymer are translations along the string such that the z-component of the center of mass is zero for all polymers, eliminating arbitrary shifts in spatial position, and translations of the overall center of mass to the coordinate origin, normalizing the position of the polymers across the simulated examples. While the first aims at recognizing the general position of the polymer with respect to the string, the latter is aimed at identifying the internal structure of globular polymers. Ge and Gi type configurations, for example, have a similar surface shape, but differ in relative position to the string and their internal crystalline structure. Finally, to help the machine recognize structural rather than size differences across all polymer types, we scaled all polymers with respect to their radius of gyration R g .

Results of unsupervised learning
In Fig. 4 we present the two-dimensional representations of the configuration space obtained from the dimensionality reduction methods mentioned above. Different columns show results after different data preprocessing steps employed (raw data without preprocesssing, standard scaling, subtracting the z-component of the center of mass, subtracting all three components of the center of mass, and data scaled by R g after subtracting the center of mass). Even without any preprocessing (leftmost column), for example, we observe that PCA can reasonably well distinguish barrel-like (B) conformations from all others. However, none of the methods can distinguish Gi, Ge, and C conformations without preprocessing. This is presumably because the z coordinates of the configurations have much higher variance than the x and y coordinates, as the system shows translational invariance in the z direction. The poor performance of ML algorithms due to different features having different scales is a common problem in machine learning. This issue can be alleviated with appropriate feature scaling techniques. Here we first test standard scaling. As the second column of Fig. 4 shows, this approach does improve the performance of the algorithms (particularly that of MDS), as a clearer separation of Gi, Ge, and C conformations can be observed. However, it is important to note that since the scaling is performed independently on individual features, these coordinate transformations lead to non-physical deformations in the polymer configurations.
A more physically intuitive scaling approach is to subtract the center of mass, which would reduce the variance of the coordinates due to the drifting of polymers in arbitrary directions. In particular, polymers have the freedom to drift along the substrate in z-direction. Therefore one would expect noticeable improvements in the results just by subtracting the z component of the center of mass alone. As the third column of Fig. 4 shows, we indeed observe improved performance in most algorithms in terms of separating previously overlapping phases observed in the analysis using raw data.
The fourth column shows the results obtained by subtracting the whole center of mass. For some algorithms (particularly MDS), subtracting all three components of the center of mass further improves phase separation. The rightmost column in Fig. 4 shows the results obtained with all coordinates furthermore normalized by scaling with the radius of gyration R g . However, we observe that Gi, Ge, and C phases are no longer distinguishable. This indicates that the length scale of polymers is a particularly important feature for distinguishing different polymer states.
In summary, we find that identifying barrel-type (B) configurations can be accomplished by all methods with suitable preprocessing steps. Telling all other structures apart is more challenging and no single scheme is able to do so alone. 2 That said, we note that the MDS method trained with data preprocessed by subtracting all three components of the center of mass seems to give the best, single overall performance, particularly since both B and C conformations are grouped into isolated clusters spatially separated from other states. Still, an overlap between Gi and Ge phases can be observed in this case since both phases differ mostly by the relative location with respect to the substrate and not in overall shape. To observe a separation between Gi and Ge structures one would need to use another procedure, like MDS or PCA with a translational normalization 2 We note that it might be possible to do so with a reduction to a 3-dimensional space though. along the z-direction only. The general finding that no one scaling approach and reduction-method combination can clearly separate all phases present in our system is probably true for other complex polymer systems as well. Depending on symmetry and specific structures, different data preprocessing and scaling methods might always have to be chosen to match all physical properties.

IV. IDENTIFYING STRUCTURAL TRANSITIONS WITH NEURAL NETWORKS
In this section we study the applicability of neural networks (NNs) to not only recognize different structures but to detect transitions points between them. While we have discussed above the desire to use large datasets for NN training in general, much more training data is potentially needed for such an endeavor since structural differences could be much more subtle between polymers close to each other in parameter space, compared to above (Sec. III) where structures are more fundamentally different from each other. To enrich our datasets we therefore apply a strategy in the spirit of oversampling augmentation [46] where we allow to record up to 100 slightly modified configurations every time the WL walker explores a low-energy region. After reaching that number, the walker has to completely "warm up" again, that is move to energies encountered well inside the random-coil phase. To ensure the data in each such batch is not effectively identical but to some degree still uncorrelated we enforce a minimum energy difference ∆E between two consecutive configurations that are added to the dataset.

A. Conventional, supervised approach
In previous research we had to rely on human intuition to define structural classes and suitable observables or order parameters to find the boundaries in parameter space between them. The structural phase diagram for low-energy states [1] (see Fig. 5 for a reduced version) was hence developed upon the ad-hoc introduction of an asymmetry parameter, for example, to locate the crossing from phase Gi to B. Such practice inevitably introduces a bias based on the human perception of structure. It is therefore, in principle, hard to judge whether or not we identified the most relevant structural features. A less biased approach that currently gets increasing attention is the use of machine learning methods to identify crossing or phase transition points between structural or thermodynamic phases [8,29,[47][48][49][50]. We here use neural networks that we train with data which can be clearly assigned to different structural classes and have them analyse polymer configurations in regions of the parameter space where such a classification is less defined.
Specifically, we investigate the transition between globular polymers absorbed to the string (Ge) and clamshell structures surrounding the string (C), two of the phases that were particularly hard to distinguish with the unsupervised methods discussed above. The NN is set up with the same parameters as above (cf. Fig. 2) with the difference that only two nodes are specified for the output layer. The network is then trained with configurations at λ = {σ f , f } = {1.0, 1.0} and λ = {2.1, 3.2}, which clearly belong to the Ge and C classes, respectively. We use the such trained network to analyse configurations at ten other parameter values in between those points (see black, diagonal line in Fig. 5) and predict their belonging to either class. When plotting the corresponding probabilities, as shown in Fig. 6 (left), we see a "crossing" of the probability curves. As one would expect, the corresponding error bars are largest around the phase intersection and decrease towards the outermost points, see Fig. 6 (right). That is, the uncertainty of the network in classifying polymer configuration is maximal around the transition from one phase to another. We assess uncertainties of the trained NN models via different methods including cross-validation [51] and queryby-committee [52]. In cross-validation, a subset of the whole dataset is held out for testing while the remaining data would be used for training. The process repeats with different held-out testing sets, resulting in a group of NN models that can be used for the estimation of statistical errors. We performed 10-fold cross-validation but it seemed to underestimate the real error of the model. It could be because the 10 resulting models are not statistically independent: each of them are trained using training datasets that overlap with each other by 80%. When these highly-correlated models are applied to make predictions on out-of-sample data (structures between λ = {1.0, 1.0} and λ = {2.1, 3.2}), they result in a small distribution (variance) around the mean, but the mean prediction might have a high discrepancy (bias) compared to the reference. Hence we report errors from query-by-committee: the whole dataset is di- vided into 10 subsets, within each subset 70% of the data were used for training and 30% of the data were used for testing. This results in 10 individually trained NN models that are truly independent and not correlated. The error bars we show in Fig. 6 therefore indicate the statistical error from multiple runs with NNs which were individually trained with independent data and also analysing different datasets. That way we capture both epistemic and aleatoric uncertainties.

B. Unsupervised: The Confusion Method
A neural network is inherently a supervised learning method and requires a dataset with preassigned labels for the adjustment of weights in the training phase. However, in certain cases it may be difficult, or even impossible, to know the correct assignment of labels beforehand. For the case of phase classification, one can circumvent this issue by identifying a window within which a given transition occurs and labeling the configurations outside this window based on the corresponding phase labels, as we did above in Sec. IV A. In particular for finite systems that do not naturally scale up to the thermodynamic limit, though, it can be challenging to reliably locate the exact point of transition this way.
The confusion method [21] provides an alternative. It not only eliminates the need for prior assignment of labels, but also results in a clearer and more precise estimate of the transition point. This method as well is a neural-network based, but semi-supervised approach for detecting phase transitions and relies on purposeful mislabeling of the data. Let λ denote a model parameter or a thermodynamic observable such as the temperature or the average energy. Assume that there exists a crit-ical point λ c at which a transition from a phase X to a phase Y occurs. When applying the confusion scheme, one first identifies a window [λ a , λ b ] within which the transition is likely to occur. Then a potential transition point λ c ∈ [λ a , λ b ] is proposed, and the label "0" (denoting phase X) is assigned to all configurations below λ c , and the label "1" (denoting phase Y) to all configurations above λ c . A neural network is then trained with this label assignment and the classification accuracy P (λ c ) obtained for a test set is recorded. This process is repeated by systematically varying λ c from λ a to λ b . The resulting curve P (λ c ) then yields a characteristic "W" shape, with the middle peak occurring at λ c = λ c [21].
This W-shaped profile of P (λ c ) can be understood as follows. For λ c = λ a , all configurations are labeled "1", and the neural network correctly predicts the assigned label for all samples, achieving 100% accuracy. Similarly, the network performs with 100% accuracy for λ c = λ b as all the configurations are labeled "0". For λ c = λ c , the assigned labels for all samples exactly match the true phase labels and, in principle, the NN can again achieve perfect accuracy. For other values of λ c , the NN sees a discrepancy between the assigned labels and the true phase labels as identified by the patterns in data. Due to this confusion, the NN learns to predict the majority label. Ultimately this yields the characteristic W shape of P (λ c ). Note that in practice this shape will likely be distorted due to finite-size effects and imperfections in the training process.
We apply the confusion method to investigate the transition from the Ge phase to the C phase along the same straight path through the phase diagram as above and indicated in Fig. 5  layer of 40 neurons. We use the same data as above and 70% of the configurations were randomly selected for training while the remaining samples were used for testing. For error estimation, the confusion scheme is repeated 10 times, each time with a different random selection of training and testing samples. Figure 7 shows the test accuracy P (λ c ) as a function of the proposed transition point λ c . The curve follows the characteristic W shape discussed above, confirming the expected transition from the Ge phase to the C phase. The middle peak indicates the location of the transition point λ c = λ c . We emphasize that the evidence for a transition provided by the confusion method is more compelling than that provided by the conventional neural network based approach discussed above in Sec. IV A. The conventional approach may falsely indicate the presence of a transition within the window [λ a , λ b ] even when all samples belong to the same phase, because also in the absence of a transition configurations may undergo slight structural changes as the parameter λ is varied between [λ a , λ b ]. More specifically, one assumes that the transition occurs within a sub-window [λ 1 , λ 2 ], and assigns the label "0" to all configurations in the interval [λ a , λ 1 ], and the label "1" to all configurations in [λ 2 , λ b ]. The neural network may then detect the gradual structural changes in the configurations as a function of λ, and establish a decision boundary between λ 1 and λ 2 such that the predictions are consistent with the assigned labels. Consequently, the curves for the averaged output neuron values may still cross each other, giving the false indication of a transition. In the confusion method however, the W-shape profile is only possible if there are abrupt, drastic changes in structure at a certain value of λ, reminiscent of a true transition. In the absence of a transition, the middle peak in P (λ c ) curve disappears, resulting in a "V" shape [21].

V. SUMMARY
In this paper we explore the applicability of various machine learning methods to recognize structures and structural transitions in a model for polymer-nanotube composites. In particular, we investigate structures that have been observed experimentally where the polymer is adsorbed at the nanotube. The two main questions we address are whether and how we can identify those structure with machine learning and how to locate the transition regions between them.
For structure recognition we test various unsupervised dimensionality reduction methods like principal component analysis or multidimensional scaling that we combine with different ways to pre-process the data. The advantage of unsupervised methods is that no prelabelling of structures is required, removing all potential human bias in structure classification. We find that while structure identification in principle is possible, no single method alone is capable of doing so. We found it particularly challenging to have the machine differentiate between globular structures where the polymer is fully wrapped around the substrate or just connects to the tube. Aside from the unsupervised methods we also employed neural network methods that do require pre-labeled input. The network was able to reliably recognize all polymer structures after suitable training.
While it is probably uncontroversial to introduce different structural phases for polymer-nanotube composites, finding the exact boundary between those phases remains a challenge since it is in general not obvious what good order parameters are. We previously introduced such parameters ad-hoc, but test here if a neural network could identify transitions between configurational phases, with and without training using configurations from the respective phases and without further human guidance or knowledge of pre-defined order parameters. This will be particularly useful since there is not a sharp, thermodynamic phase transition. We find that neural-network methods still indicate a transition, most notably the confusion method. However, since such structural transitions of finite systems potentially happen in different steps and over a broader region in parameter space, different machine learning methods or neural networks might pick up different steps in this transition at slightly different parameter values. In that sense, results for the crossing point shown in Figs. 6 and 7 are not necessarily contradicting, when keeping also in mind that data has to be binned for the confusion method, leaving a corresponding uncertainty in the exact position of the crossing. That said, we also note that the traditional method of training the network with labeled configurations from both phases has to be used with care since it can potentially detect a crossing even if there was no phase transition. The main advan-tage of the confusion method here is therefore that it provides evidence for a transition between Ge and C. Otherwise, the shape of the detection accuracy graph would be a "V"-shape rather than a "W"-shape.
Overall, we confirm that defining structural transitions in our system is reasonable, in principle. We also conclude, though, that we might not have been successful initially [1] in finding the best order parameter for all transitions, in particular for the crossing between "Ge" and "C" structures, as evidenced by the results shown in Figs. 5-7. In that sense, the machine learning methods applied here can be a valuable complement to more conventional methods of detecting structural transitions used earlier, as they remove the necessity of identifying or defining explicit order parameters beforehand and therefore provide a potentially less biased approach to structure recognition and classification.


Training Diverse High-Dimensional Controllers by Scaling Covariance Matrix Adaptation MAP-Annealing

Pre-training a diverse set of robot controllers in simulation has enabled robots to adapt online to damage in robot locomotion tasks. However, finding diverse, high-performing controllers requires specialized hardware and extensive tuning of a large number of hyperparameters. On the other hand, the Covariance Matrix Adaptation MAP-Annealing algorithm, an evolution strategies (ES)-based quality diversity algorithm, does not have these limitations and has been shown to achieve state-of-the-art performance in standard benchmark domains. However, CMA-MAE cannot scale to modern neural network controllers due to its quadratic complexity. We leverage efficient approximation methods in ES to propose three new CMA-MAE variants that scale to very high dimensions. Our experiments show that the variants outperform ES-based baselines in benchmark robotic locomotion tasks, while being comparable with state-of-the-art deep reinforcement learning-based quality diversity algorithms. Source code and videos are available at https://scalingcmamae.github.io

I. INTRODUCTION
By generating a diverse collection of controllers, we can endow a robot with a wide variety of useful behaviors, such as moving an arm into different configurations [1], [2], walking in every direction [3], [4], pushing an object through many trajectories [5], or traveling to any position in a maze [6], [7], [8]. Having a collection of controllers has proven particularly useful in robotic locomotion, where the collection enables a walking robot to adapt to damage [9], [10], [11]. The controllers differ by how long each foot contacts the ground, such that if one foot is damaged, the robot can select a controller that does not rely on that foot.
Searching for these diverse controllers may be viewed as a quality diversity (QD) optimization problem [6]. In QD, the goal is to find solutions φ that are diverse with respect to one or more measure functions m(φ) while maximizing an objective function f (φ). In robotic locomotion, we search for controller policies π φ parameterized by φ. Each controller should satisfy a unique output of the measure function by using its feet in a different manner from the other controllers, while optimizing the objective by walking forward quickly.
To solve a QD problem, a QD algorithm must balance two aspects given a limited compute budget: exploring the measure space and optimizing the objective. Furthermore, one of these aspects may be more important for a given problem. For example, in locomotion, exploring corresponds to finding new controllers which use the robot's feet a different amount, while optimizing corresponds to making existing controllers walk faster. Prior work [12] shows that Fig. 1: We propose variants of the CMA-MAE algorithm which scale to high-dimensional controllers. The variants maintain a Gaussian search distribution with mean φ * and approximate covariance matrixΣ. Solutions φ i sampled from the Gaussian are evaluated and inserted into an archive, where they generate improvement feedback ∆ i based on their objective value f (φ i ) and a threshold t e which each archive cell maintains. Finally, the Gaussian is updated with an evolution strategy (ES). Our variants differ from CMA-MAE by selecting scalable ESs, as the CMA-ES used in CMA-MAE has Θ(n 2 ) time complexity per sampled solution.
optimization is more difficult than exploration in this domain, as it is trivial to explore the measure space by standing in place and lifting the feet up and down to achieve a new style of foot usage.
Prior algorithms [13], [12] seem to strike a balance between these two aspects of QD in locomotion tasks, leading to state-of-the-art results. However, these algorithms have practical limitations due to their dependence on deep reinforcement learning (RL) methods. Namely, they require a GPU, must perform time-consuming training of a neural network, and have a large number of hyperparameters.
Recent work [14] suggests evolution strategies (ES) as a compelling alternative to deep RL methods when optimizing a single controller. Compared with deep RL, ESs can run entirely on CPU and do not require network training, and ESs such as CMA-ES [15] are designed to have almost no hyperparameters. Given these benefits, prior work [12], [11] has developed QD algorithms based on ESs, but these methods have not yet been able to match the performance of deep RL-based QD methods.
On the other hand, the recently proposed ES-based Covariance Matrix Adaptation MAP-Annealing (CMA-MAE) algorithm [16] has proven adept at trading off between the exploration and objective optimization aspects of QD. Through a single hyperparameter α ∈ [0, 1], CMA-MAE can be tuned to optimize only the objective (as in CMA-ES), or prioritize moving away from previously discovered high-quality solutions (as in CMA-ME [17]), and anywhere in between. Selecting any value of α strictly between 0 and 1 yields theoretical benefits that result in state-of-theart performance on several QD benchmark domains. We hypothesize that because CMA-MAE balances this tradeoff, the algorithm would excel in robotic locomotion tasks with hard-to-optimize objectives.
While CMA-MAE excels at moderate dimensional domains, the algorithm is computationally intractable for robotic locomotion due to the high dimensionality of modern neural network controllers. CMA-MAE relies on CMA-ES [15], whose Θ(n 2 ) time complexity per sampled solution, where n is the number of parameters, makes it impossible for CMA-MAE to optimize the thousands or even millions of parameters found in such controllers.
CMA-ES's complexity arises from how it models the search distribution by maintaining a Gaussian with a full rank n×n covariance matrix. However, a full rank covariance matrix is often not necessary, so prior work [18], [19], [20], [21], [22], [23] scales CMA-ES to higher dimensions by replacing the covariance matrix with sparse approximations. These CMA-ES variants scale to solve a variety of highdimensional optimization problems.
Thus, our key insight is that we can scale CMA-MAE to high-dimensional controllers by adopting the same approximations in its CMA-ES components. Though our CMA-MAE variants leverage sparse covariance matrix approximations, we hypothesize these variants will be competitive.
We scale CMA-MAE by proposing three CMA-MAE variants, which we evaluate on robotic locomotion tasks from QDGym [24]. We show that our variants are the highestperforming QD methods based solely on ES. Furthermore, they inherit the aforementioned practical benefits of ES, and they are comparable to the state-of-the-art deep RL-based QD method PGA-ME [13] on three of four tasks. We are excited about future applications in other domains, such as robotic manipulation [25] and scenario generation [26].
II. PROBLEM STATEMENT Quality diversity (QD). Drawing from recent work [27], we define QD as considering an objective function f (φ) and k-dimensional measure function m(φ), 1 where φ ∈ R n is an n-dimensional solution. The outputs of m form a kdimensional measure space X . The QD objective is to find solutions φ such that ∀x ∈ X , m(φ) = x and f (φ) is maximized. Solving this QD objective would require infinite memory because X is a continuous space, so algorithms based on MAP-Elites [9], [10] relax the QD objective by discretizing X into a tesselation Y of M cells. Then, the QD objective is to find an archive (i.e., a collection) A containing M elites φ i , where each φ i has measures m(φ i ) corresponding to a cell in Y and f (φ i ) is maximized. We express this relaxed QD objective as max Quality diversity for reinforcement learning (QD-RL). As defined in prior work [12], QD-RL is a special instance of QD in which φ parameterizes a reinforcement learning (RL) agent's policy π φ and the objective is the expected discounted return of the agent. QD-RL extends Markov Decision Processes (MDPs) [28] and is formulated as a tuple (S, U, p, r, γ, m). Here, on each timestep of an episode, an agent observes state s ∈ S and takes action a ∈ U with policy π φ (a|s). The agent transitions to the next state s ∈ S according to s ∼ p(·|s, a) and receives reward r(s, a, s ). The expected discounted return is defined as f (φ) = E T t=0 γ t r(s t , a t , s t ) for a discount factor γ ∈ (0, 1) and episode length T . Finally, QD-RL includes a kdimensional measure function m(φ) which describes the agent's behavior during an episode.

III. BACKGROUND A. MAP-Elites
Many QD algorithms, including the ones in this paper, are based on the MAP-Elites algorithm. The original MAP-Elites algorithm [9], [10] divides the measure space equally along each dimension, resulting in an archive consisting of grid cells. Initially, MAP-Elites samples solutions from a predefined distribution and inserts them into the archive by mapping each solution to a cell based on its measures. On subsequent iterations, MAP-Elites samples a solution φ from the archive, mutates it by adding Gaussian noise ∼ N (0, I), and adds the new solution φ to the archive. Crucially, if a solution ever lands in the same cell as a previous solution, it is only stored in the archive if it has a higher objective. Hence, solutions in the archive are referred to as elites.

B. Large-Scale Evolution Strategies (ES)
ESs [41] optimize continuous parameters by adapting a population of solutions, such that the population is more likely to have high performance. Here we discuss ESs designed for high-dimensional search spaces. In Sec. IV, we show how these ESs help scale CMA-MAE.
Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [15]. CMA-ES is an approximate second-order method which follows a natural gradient [42] and demonstrates state-of-the-art results in black-box optimization [43].
CMA-ES maintains a Gaussian N (µ, Σ) that models the distribution of search directions. In every iteration, CMA-ES samples λ solutions from the Gaussian and updates it based on rankings of the solutions' performance.
CMA-ES requires Θ(n 2 ) space and has Θ(n 2 ) runtime per sampled solution due to its n × n covariance matrix Σ. The runtime stems from two operations. First, updating Σ requires matrix-vector multiplications. Second, since it is easy to sample from the standard Gaussian N (0, I) on a computer, sampling from the Gaussian is implemented as Computing the transformation matrix Σ iterations. Nevertheless, this complexity is intractable for high n. To scale CMA-ES, prior works [18], [19], [20], [21], [22], [23] assume that a full Gaussian is unnecessary for modeling the search distribution in high dimensions and thus replace the covariance matrix with a more efficient approximation. We discuss two such approaches below.
Separable CMA-ES (sep-CMA-ES) [18]. As the simplest large-scale CMA-ES variant, sep-CMA-ES assumes that the problem is separable, i.e., that it may be solved by optimizing each parameter individually. Hence, sep-CMA-ES constrains the covariance matrix to be diagonal, so it only stores n parameters. Updates and sampling then require Θ(n) runtime per solution, as the matrix-vector multiplications are replaced with dot products and Σ 1 2 in Eq. 1 is computed by taking the square root of the diagonal elements. sep-CMA-ES outperforms CMA-ES on separable benchmarks and even on some non-separable problems.
Limited-Memory Matrix Adaptation ES (LM-MA-ES) [19]. LM-MA-ES is another large-scale CMA-ES variant. It maintains a rank-k approximation of the transformation matrix Σ 1 2 . This approximation consists of k vectors which represent different directions of the search distribution. Since the approximation is rank-k, the space complexity is Θ(kn), and the time complexity is Θ(kn) per solution. Overall, LM-MA-ES performs well on optimization benchmarks, matching other large-scale versions of CMA-ES.
OpenAI-ES [14]. Unlike the algorithms above, OpenAI-ES is not a CMA-ES variant. Instead, inspired by Natural Evolution Strategies [44], OpenAI-ES maintains an isotropic Gaussian with constant covariance matrix σI, and it only updates the mean φ * via gradient ascent with Adam [45]. The gradient for this update is computed based on the ranknormalized [46], [44] performances of λ solutions sampled from the Gaussian with mirror sampling [47]. OpenAI-ES has performed well across RL benchmarks [14], [48], [49]. It has a time complexity of Θ(n) per sampled solution.

C. Combining ES with MAP-Elites
Here we discuss algorithms which enhance MAP-Elites by integrating ESs.
MAP-Elites with Evolution Strategies (ME-ES) [11]. ME-ES maintains a solution point φ * and modifies it with OpenAI-ES in two phases: exploit, where OpenAI-ES optimizes φ * to have high objective value, and explore, where OpenAI-ES optimizes φ * to be novel [7], [50] (far away from previous solutions in measure space). Both phases begin by sampling a new φ * from the archive. They then run for k iterations, with only φ * inserted into the archive after each iteration. ME-ES is a baseline in our experiments.
Covariance Matrix Adaptation MAP-Annealing (CMA-MAE) [16]. One promising class of QD algorithms combines CMA-ES [15] with MAP-Elites. We directly extend the most recent such algorithm, CMA-MAE, which has shown state-of-the-art performance in QD benchmark domains. CMA-MAE maintains an archive which records a threshold t e for each cell e. t e is initialized to a user-defined minimum objective min f . When a solution φ is inserted into the archive, it is placed into its corresponding cell e if its objective value f (φ) exceeds t e . Then, t e is updated via polyak averaging, i.e. t e ← (1 − α)t e + αf (φ), where α ∈ [0, 1] is the archive learning rate. Note that the solution's objective value only needs to cross the threshold, rather than the objective value of the elite previously in the cell. Hence, the archive does not always store the best solutions like MAP-Elites does, so implementations must track these best solutions separately.
For generating solutions, CMA-MAE queries a set of emitters on each iteration. Each emitter maintains its own Gaussian search distribution with mean φ * and covariance matrix Σ. In every iteration, the emitter samples solutions from this Gaussian and inserts the solutions into the archive. When a solution φ is inserted, the archive returns an improvement value ∆ i ← f (φ) − t e , i.e., the difference between the objective value of the solution and its corresponding cell's threshold. The emitter ranks its generated solutions by these improvement values. Based on this ranking, the emitter updates φ * and Σ with CMA-ES, such that it is more likely to generate solutions with high improvement in the future.
A key parameter in CMA-MAE is the archive learning rate α. When α = 0, the threshold always remains at min f , and hence the improvement is always equivalent to the objective value (minus the constant min f ). This makes CMA-MAE equivalent to CMA-ES, as it optimizes solely for the objective function. However, when α = 1, the threshold is always equivalent to the objective value of the elite currently in the cell, which means there is minimal improvement for inserting a solution into that cell. In this case, CMA-MAE is equivalent to the previously developed CMA-ME [17], where discovering new cells in the archive (i.e., exploration) always takes priority over improving the objective value of cells where there is already a solution. Hence, by varying α from 0 to 1, we trade off between optimizing the objective and moving away from previously discovered elites. In particular, we use a low α value due to the need for optimization in robotic locomotion domains.
Covariance Matrix Adaptation MAP-Elites via a Gradient Arborescence (CMA-MEGA) [27]. In the differentiable quality diversity (DQD) setting, where exact objective and measure gradients are available, CMA-ME and gradient ascent have been combined to form the state-of-the-art CMA-MEGA algorithm. To extend CMA-MEGA to QD-RL, prior work replaces the exact gradients which CMA-MEGA relies on with gradient approximations, since the RL environment is usually non-differentiable. This results in two variants: CMA-MEGA (ES) approximates all gradients with OpenAI-ES [14], while CMA-MEGA (TD3, ES) approximates the objective gradient with TD3 [40] and measure gradients with OpenAI-ES. These variants are baselines in our experiments.

IV. SCALING COVARIANCE MATRIX ADAPTATION MAP-ANNEALING
In CMA-MAE, each emitter relies on CMA-ES to update the Gaussian search distribution from which it samples solutions. Since CMA-ES requires Θ(n 2 ) runtime per solution (where n is the solution dimension), CMA-MAE is intractable for neural network controllers, where there are thousands to millions of parameters. To scale CMA-MAE, we propose three variants which substitute CMA-ES with large-scale ESs. These variants differ primarily in the complexity of their covariance matrix approximation: • sep-CMA-MAE substitutes sep-CMA-ES [18], which constrains the covariance matrix Σ to be diagonal, resulting in Θ(n) complexity. • LM-MA-MAE substitutes LM-MA-ES [19], which approximates the transformation matrix Σ 1 2 with k n n-dimensional vectors, leading to Θ(kn) complexity. • OpenAI-MAE substitutes OpenAI-ES [14], which sets the covariance matrix Σ to be the identity matrix. Though Σ is constant, various operations still necessitate Θ(n) complexity. The listed complexities refer to (1) the space required per emitter, as each emitter maintains its own ES instance, and (2) the runtime required per sampled solution, which is the same regardless of the number of emitters.
Algorithm 1 shows the pseudocode for these CMA-MAE variants, and Fig. 1 shows an overview of the variants. Each variant proceeds as follows. First, it initializes the archive along with each emitter's ES parameters (line 2-3). This step includes initializing the covariance matrix approximationΣ in lieu of CMA-MAE's full covariance matrix Σ. Next, the variant repeatedly queries the emitters for solutions (line 4). Each emitter (line 5) samples λ solutions from the distribution N (φ * ,Σ) (line 6-7). Note that the sampling procedure depends on the approximation employed by the variant. Once sampled, the solutions are evaluated and inserted into the archive if they cross their cell's threshold t e (line [8][9][10][11][12][13][14][15]. Then, the ES's parameters (includingΣ) are updated based on the improvement ranking of the solutions (line [16][17], such that the emitter is more likely to sample solutions with high improvement on the next iteration. Finally, the emitter is restarted if the ES converges (line [18][19][20]. We adopt the default update and convergence rules from each ES.

V. EXPERIMENTS
To test our CMA-MAE variants' ability to create diverse, high-performing locomotion controllers, we evaluate it on Replace the solution in cell e of archive A with φ i four robotic locomotion tasks from the QDGym [24] benchmark based on PyBullet [51], [52] and OpenAI Gym [53]. In each environment (Fig. 2), the objective is to walk forward quickly, while the measures track the proportion of time in an episode that each of the robot's feet touches the ground, i.e., if a robot has four legs, it has four measures.
As prior work [12] notes, these environments are challenging because they require focusing on objective optimization. It is easy to fill the archive with poor controllers that achieve different measures by standing still and lifting the feet up and down, but making these controllers walk forward is much more difficult. We compare our variants with five baselines: PGA-ME, two CMA-MEGA variants (CMA-MEGA (ES) QD Ant QD Half-Cheetah QD Hopper QD Walker Fig. 2: QDGym locomotion environments [24]. Images reproduced with permission from prior work [12].

A. Hyperparameters
As in prior work [13], [12], the robot controller is a neural network mapping states to actions. The network has two hidden layers of size 128 and tanh activations on every layer but the input, and it is initialized with Xavier initialization. The number of parameters in each domain is as follows: QD Ant (21,256), QD Half-Cheetah (20,742), QD Hopper (18,947), QD Walker (20,230).
All algorithms are allocated 1 million evaluations, with each solution evaluated for one episode. PGA-ME, CMA-MEGA (ES), CMA-MEGA (TD3, ES), and ME-ES use the same hyperparameters as in their prior works [13], [12], [11], except ME-ES uses 200 samples. MAP-Elites uses a noise standard deviation of σ = 0.02 and a batch size of 100.
All CMA-MAE variants run with a low archive learning rate α = 0.001 to encourage optimization of the objective and with ψ = 5 emitters, each with batch size λ = 40 and initial step size σ = 0.02. LM-MA-MAE sets k = λ = 40. The Adam optimizer for OpenAI-ES in OpenAI-MAE uses learning rate 0.01 and L2 regularization coefficient 0.005.

B. Experimental Design
We conduct a between-groups study with the algorithm (LM-MA-MAE, sep-CMA-MAE, OpenAI-MAE, PGA-ME, CMA-MEGA (ES), CMA-MEGA (TD3, ES), ME-ES, MAP-Elites) and environment (QD Ant, QD Half-Cheetah, QD Hopper, QD Walker) as independent variables and QD score as dependent variable. The QD score [6] measures performance holistically by summing the objectives of all elites in the archive. To ensure that no elite subtracts from the QD score (objectives may be negative), we subtract min f from all elites' objectives before computing the QD score. We repeat our experiments 5 times to test three hypotheses: H1: All CMA-MAE variants will outperform prior ESbased methods, i.e., ME-ES and CMA-MEGA (ES).

H3:
The performances of the CMA-MAE variants will be ordered as LM-MA-MAE > sep-CMA-MAE > OpenAI-MAE.
H1 and H2 are based on the key difference between the variants and the baselines, i.e., the smooth improvement ranking based on archive thresholds. Since CMA-MEGA (ES) and CMA-MEGA (TD3, ES) use a standard MAP-Elites archive (equivalent to setting α = 1 in CMA-MAE), we think they will focus too much on exploration and be unable to optimize existing solutions. Meanwhile, PGA-ME and ME-ES separate measure space exploration from objective optimization; PGA-ME has two types of operations, while ME-ES has two distinct phases. We hypothesize that this separation will not work as well as blending the two aspects together via the improvement ranking.
The ordering in H3 is based on the complexity of each variant's covariance matrix approximation. Essentially, given the difficulty of the locomotion domains, we hypothesize that a better representation of the search distribution will lead to higher performance.

C. Implementation
Each algorithm runs with 100 CPUs on a highperformance cluster. Algorithms which train TD3 are allocated one NVIDIA Tesla P100 GPU. We obtained data for all baselines with consent from the authors of prior work [12]. All algorithms are implemented with the pyribs library [55].

VII. DISCUSSION AND CONCLUSION
Our results show that ES-based methods are a viable alternative in quality diversity for reinforcement learning (QD-RL). First, the analysis of H1 shows that the CMA-MAE variants are the highest-performing ES-based method in QD-RL, either matching or exceeding the existing CMA-MEGA (ES) and ME-ES on all domains. Additionally, the analysis of H2 shows that the variants are comparable to RL-based methods (including the state-of-the-art PGA-ME) in 3 of 4 tasks. While the variants do not solve QD Walker, prior work [12] highlights the importance of deep RL in this task, as only algorithms with TD3 (PGA-ME and CMA-MEGA (TD3, ES)) have performed well here. Overall, these results show that the improvement ranking and annealing mechanisms of CMA-MAE enable it to focus on the optimization of the objective, which is necessary in the robotic locomotion tasks.
We note that a benefit of ES-based methods is their generality. While PGA-ME performs well on the robotic locomotion tasks, it is designed only for QD-RL settings, since its TD3 component is applicable only in MDP domains. On the other hand, the CMA-MAE variants treat the objective as a black box and can be applied to general QD problems outside of robotics.
Furthermore, ES-based methods have attractive practical benefits. Specifically: (1) The CMA-MAE variants are noticeably light on computation. While PGA-ME and CMA-MEGA (TD3, ES) both train a TD3 critic, a lengthy process which requires a GPU, the CMA-MAE variants run entirely on CPU and are three times faster (7 hrs vs 20 hrs in QD Ant and QD Half-Cheetah, and 4 hrs vs 13 hrs in QD Hopper and QD Walker).
(2) The CMA-MAE variants have substantially fewer hy-perparameters than PGA-ME and the CMA-MEGA variants since they depend on CMA-ES and its variants, which are designed to be parameterized only by an initial step size σ and batch size λ. Hence, the CMA-MAE variants only require 5 hyperparameters (ψ, λ, σ, α, min f , see Algorithm 1). In contrast, PGA-ME and CMA-MEGA (TD3, ES) require a large number of parameters for TD3 and their other components.
When comparing the CMA-MAE variants in H3, we observe that OpenAI-MAE performs significantly worse than one other variant in three out of four tasks. This is because OpenAI-MAE approximates the covariance matrix with a constant identity matrix, as opposed to the more complex approximations of sep-CMA-MAE and LM-MA-MAE. On the other hand, there does not seem to be a significant difference between sep-CMA-MAE and LM-MA-MAE. This result warrants further investigation, and it indicates that both the rank-k and diagonal approximations of the covariance matrix are sufficient for modeling the distribution of search directions in the locomotion tasks.
In conclusion, our work scales CMA-MAE to robotic locomotion controllers by replacing the CMA-ES component of CMA-MAE with efficient approximations. Our results show the CMA-MAE variants are effective in benchmark locomotion tasks and have fewer hyperparameters than the deep RL baselines, making them attractive to practitioners. Since CMA-MAE is a general purpose algorithm, we are excited about future work that will test the CMA-MAE variants in domains beyond robot locomotion, such as robotic manipulation [25] and scenario generation [26]. Future work will also test the pretrained controllers in the real world and will explore the computational benefits of recent hardwareaccelerated frameworks [56].


chowdsp_wdf: An Advanced C++ Library for Wave Digital Circuit Modelling

chowdsp_wdf is a C++ library for implementing real-time wave digital models of analog circuits. chowdsp_wdf differs from existing wave digital modelling libraries by providing a template meta-programming interface for modelling circuits with a fixed topology, and providing support for explicit SIMD acceleration. The motivation and design of the library are described, as well as real-world use-cases, and performance comparisons with other wave digital modelling libraries.

INTRODUCTION
Wave Digital Filters (WDFs) are a circuit modelling paradigm originally developed by Alfred Fettweis in the 1970's and '80's [1]. WDFs have gained popularity in recent years for developing "virtual analog" (VA) emulations of audio circuits [2]. VA emulations are often implemented as part of software synthesizers or audio effects that are designed to run in real-time, so the run-time performance of these emulations is an important consideration.

Wave Digital Filters Background
While this paper will not attempt to provide a complete introduction to the WDF paradigm, there are several properties of WDFs which should be mentioned. First, WDFs are "modular," meaning that a model of a circuit can be separated into its component parts. This property makes it possible for a wave digital modelling library to implement a single "Resistor" element, which can then be re-used in any wave digital circuit model (and so on with most other circuit elements). These element models may also provide interfaces for changing their parameters (e.g. the resistance value of a resistor).
Wave digital models are comprised of the circuit's element models, along with a set of "adaptors" describing how the circuit elements are connected. Common adaptors include "series" and "parallel" adaptors, as well as polarity inverters. The circuit elements are typically arranged in a tree-like structure, so that each adaptor has a single "parent" element and one or more "child" elements. When processing signals through the wave digital model, each element requires only information that may be provided by it's child elements. However, when the circuit topology changes, or when a component value is changed (e.g. a potentiometer has been adjusted), it is sometimes necessary for an element to signal the change to its parent element.
Finally, rather than using traditional circuit quantities such as voltage and current, WDFs operate on wave variables. Most WDF circuit elements accept a single "incident" wave and output a single "reflected" wave. Adaptors accept N incident waves and output N reflected waves, where N is the total number of child and parent elements connected to the adaptor. The process of computing the output reflected waves from the incoming incident waves for an adaptor may generally be expressed in the form of a scattering matrix, however one-multiply forms of the scattering operation exist for series and parallel adaptors.

PREVIOUS WORK
Several libraries for creating real-time wave digital models exist, including RT-WDF, a C++ library published in 2016 [3], and wdmodels, a Faust library published in 2021 [4]. There are also several libraries not being considered in this paper (e.g. [5]) since they are written in languages which are not well-suited for real-time audio applications, however, these libraries may be useful for prototyping and for performing non-real-time simulations.
RT-WDF is based on C++'s concept of run-time polymorphism. Each wave digital element and adaptor is derived from a base class which provides interfaces for all shared WDF functionality. Additionally, RT-WDF defines the sample rate as a parameter of the constructor for classes that require knowledge of the system sample rate (capacitors, inductors, etc.). Since the system sample rate may not be known at compile-time, circuit elements are typically constructed dynamically using "heap" memory, i.e. with C++'s std::unique_ptr<>. For custom adaptors that require more complex scattering matrices, RT-WDF interfaces with the Armadillo 1 library to perform optimized matrix operations.
wdmodels is written using the Faust programming language. Faust is a functional programming language for expressing audio processing algorithms, which may be compiled into C++, LLVM Bitcode, Rust, and other programming languages. Faust requires the signal flow of the system to be fixed at compilation time, meaning that the topology of the circuit being modelled must be fixed at compile-time as well. Faust also does not provide the programmer with mechanisms for interacting with hardware-level instructions, including SIMD intrinsics (helpful for implementing faster matrix operations for adaptors that require more complex scattering matrices), or bit-twiddling (helpful for implementing approximations of various math functions, see e.g. [6]).

LIBRARY DESIGN
chowdsp_wdf intends to provide the best features of both RT-WDF and wdmodels, with improved performance and flexibility. 2 The library is open-source and is published under the BSD "3-clause" license.

Run-Time Polymorphism
For situations where the circuit topology may need to be changed at run-time, chowdsp_wdf adopts a similar strategy as RT-WDF, using C++'s run-time polymorphism, except that components may have their sample rate changed after construction. This property allows circuit elements to be constructed using "stack" memory rather than "heap" memory, resulting in improved locality of reference [7]. However, the polymorphic approach has significant performance limitations when compared with wdmodels's strategy of defining the circuit topology at compile-time (limiting though it may be). Whenever an adaptor in RT-WDF needs to access some property of one of its children, it must do so through the polymorphic interface of the child's base class, which is typically implemented as a virtual method table or "vtable" [8]. Performing an operation by vtable lookup introduces a small performance overhead and prevents the compiler from performing an inline expansion of the code being executed [9].
The negative performance effects of implementing a wave digital circuit model using run-time polymorphism can be demonstrated by analyzing the assembly instructions generated by an optimizing compiler. When compiling the simple circuit model shown in listing 1 using the Clang 13.1.6 compiler, with the -std=c++20 -O3 compiler flags, the generated assembly for the process() method contains four virtual function call instructions (two for S1 to receive the reflected waves from R1 and R2, and two to send incident waves to those same elements). These virtual function calls add a small amount of computational overhead and prevent the compiler from performing an inline expansion of the code being run behind the call instructions, which then prevents further optimizations from taking place. For larger circuit models containing more circuit elements, these issues become an increasingly significant performance bottleneck.

Optimizations Via Template Meta-Programming
chowdsp_wdf provides an API for constructing wave digital circuit models using compile-time template meta-programming. While this API does not have the flexibility to construct arbitrary circuit models at run-time, it can provide significant performance improvements for cases where the circuit topology is fixed at compiletime. The fundamental idea behind this design choice is that giving the compiler the maximum possible information about the circuit model will allow the compiler to perform the best possible optimizations for that model. In the compile-time API, resistors, capacitors, and most other "one-port" circuit elements are defined almost identically to their run-time counterparts. However, compiletime adaptors are defined with additional template arguments to determine the types of the child elements which are to be connected to the adaptor. Since the adaptor is aware of the types of its child elements at compile-time, all audio processing operations can be performed directly, i.e. without requiring a vtable lookup. Listing 2 shows the same circuit model as listing 1, implemented using chowdsp_wdf's compile-time API. When compiling this code with the same compiler and settings, the generated assembly contains zero call instructions, implying that all the necessary interfaces between the adaptor and its child elements have been inlined by the compiler. A further comparison between the performance of chowdsp_wdf's run-time and compile-time API's is provided in Section 4.

Data Type Abstraction
RT-WDF uses double-precision floating point numbers to store all quantities in the wave digital circuit model. As with all Faust code, wave digital models written with wdmodels use singleprecision floating point values by default, however, the Faust compiler contains optional arguments for using double-precision or quad-precision floating point values instead.
chowdsp_wdf provides a template abstraction for the data type used to store quantities in the circuit. By default, chowdsp_wdf supports C++'s native floating point data types, as well as SIMD wrappers around those floating point data types, via the XSIMD C++ library. 3 The ability to construct wave digital circuit models using SIMD data types creates several interesting optimization opportunities for implementers of wave digital circuit models.

Polyphonic Synthesizer
Consider a polyphonic synthesizer being implemented for a platform that supports Intel SSE or ARM NEON SIMD intrinsics. The synthesizer voices could be implemented using a SIMD data type containing 4 single-precision floating point numbers, providing up to a 4x performance improvement over the same synthesizer using scalar floating point numbers for each voice. Listing 3 shows a minimal example of how such a synthesizer voice might be implemented.

Parallel Circuits
Constructing wave digital circuit models using SIMD data types may also be useful when emulating devices that contain multiple instances of the same sub-circuit. One example of this phenomenon is the Buchla 259 wavefolder [10], which contains 5 instances of the same "folder cell" circuit in parallel with each other. A similar approach could also be used to emulate bucket-brigade device circuits, which contain many capacitor "bucket" circuits in series. For example, a bucket-brigade device with 1024 buckets could be split into 4 parallel sub-circuits each containing 256 buckets, arranged such that the output of the first sub-circuit feeds into the input of the second sub-circuit and so on.

R-Type Adaptors
For circuits with topologies that cannot be broken down strictly into series and parallel connections between components, constructing a wave digital model of the circuit may require the use of an R-Type adaptor [11]. A simple R-Type adaptor may have any number of child elements, and up to one parent element, and uses a scattering matrix to compute the outgoing reflected waves from the incoming incident waves.
chowdsp_wdf contains interfaces for constructing and working with R-Type adaptors, based on the original implementation provided in [12]. In order to update the scattering matrix when a "downstream" component value changes, the R-Type adaptors in the compile-time API are implemented using a "template functor" so that the adaptor can use some custom logic to compute the scattering matrix. Listing 4 demonstrates a simple circuit model containing an R-Type adaptor with 3 child elements and no parent elements. For performing operations with the scattering matrix chowdsp_wdf will optionally use XSIMD to perform SIMDaccelerated matrix operations. Listing 4: A simple circuit model containing an R-Type adaptor.

Deferring Adaptor Updates
For circuits with control parameters, a situation may arises when several component values need to be updated at once. In this situation, some computations may become redundant if an adaptor is alerted multiple times of some changes in "downstream" circuit elements. This redundancy can affect run-time performance, particularly for adaptors that require relatively large computations in response to changes in downstream elements, as is often the case with R-Type adaptors. To mitigate these performance impacts, chowdsp_wdf provides a mechanism for "deferring" adaptor updates until after all downstream elements have been updated.

PERFORMANCE COMPARISONS
In order to compare the run-time performance of each modelling library, a variety of circuit models were implemented with each library. 4 Chosen circuits include a passive second-order RC lowpass filter (LPF2), a sub-circuit from the Klon Centaur guitar pedal (FF-2), an RC diode clipper, and the tone stack circuit from the '59 Fender Bassman guitar amplifier [13]. The LPF2 and FF-2 circuits were chosen as example circuits of varying complexity with wave digital models containing only linear elements and simple adaptors. The diode clipper circuit was chosen as an example circuit containing a single nonlinear element. The Bassman tone stack was chosen as an example circuit requiring an R-Type adaptor. No implementation of the diode clipper circuit was attempted with RT-WDF, since the library does not provide wave domain implementations of diode elements. Each model was also implemented using chowdsp_wdf's compile-time API, as well as chowdsp_wdf's polymorphic run-time API.
A 1000 second long audio signal at a sample rate of 48 kHz was processed through each model implementation, using a 2018 Mac Mini, with a 3.2 GHz Intel Core i7 CPU. The time needed for each model implementation to process the signal was measured, and is shown in table 1. For each circuit model, the chowdsp_wdf polymorphic implementation outperforms the RT-WDF implementation, with the performance difference widening for the more complicated circuits. Between the two model implementations that are fixed at compile-time, chowdsp_wdf's compile-time API outperforms the wdmodels implementation for all circuit models, excepting the FF-2 circuit.
The largest difference in performance between the chowdsp_wdf compile-time API and the wdmodels library can be seen in the diode clipper circuit. Both libraries use diode implementations based on the Werner et al.'s model [14], which requires evaluation of the Lambert W function. The large performance difference can likely be ascribed to the fact that chowdsp_wdf implements the Lambert W with a numerical approximation that utilises floatingpoint bit twiddling [6], while wdmodels uses a Newton-Raphson solver with a fixed number of iterations.

EXAMPLES
For students and engineers looking to learn how to implement circuit models using chowdsp_wdf, a GitHub repository 5 has been created containing models of several simple circuits (voltage divider, current divider, RC lowpass filter, diode clipper, etc.), as well as several more complex audio circuits, such as the Baxandall tone control circuit [15], and the Sallen-Key filter [16]. There are also several open-source examples of chowdsp_wdf being used in audio plugins. 6 ,7

CONCLUSION
chowdsp_wdf is a C++ library for implementing wave digital circuit models, with separate APIs for construct circuit models with run-time polymorphism, or with compile-time template metaprogramming. The library also supports explicit SIMD acceleration, R-Type adaptors, and deferred adaptor updates. Circuit models implemented with chowdsp_wdf typically exhibit superior run-time performance when compared with the same circuit model implemented using other wave digital modelling libraries.
Future work will involve setting up bindings to use the chowdsp_wdf library from other programming languages, including Rust, and Javascript (through the web audio API). Additionally, the author will continue to investigate further performance improvements for the library as whole.

ACKNOWLEDGMENTS
The author would like to thank Jingjie Zhang, Kurt Werner, and Dirk Roosenburg for sharing their knowledge about the theory and implementation of wave digital filters, as well as Sam Schachter for his help with implementing R-Type adaptors. Thanks as well to Eyal Amir and Paul Walker for their crucial insights regarding C++ optimization.


extfieldsofstudy:
Computer Science
Engineering

text:
Analyzing Acoustic Word Embeddings from Pre-trained Self-supervised Speech Models

Given the strong results of self-supervised models on various tasks, there have been surprisingly few studies exploring self-supervised representations for acoustic word embeddings (AWE), fixed-dimensional vectors representing variable-length spoken word segments. In this work, we study several pre-trained models and pooling methods for constructing AWEs with self-supervised representations. Owing to the contextualized nature of self-supervised representations, we hypothesize that simple pooling methods, such as averaging, might already be useful for constructing AWEs. When evaluating on a standard word discrimination task, we find that HuBERT representations with mean-pooling rival the state of the art on English AWEs. More surprisingly, despite being trained only on English, HuBERT representations evaluated on Xitsonga, Mandarin, and French consistently outperform the multilingual model XLSR-53 (as well as Wav2Vec 2.0 trained on English).

INTRODUCTION
Speech tasks such as query-by-example, voice search, keyword spotting, and word discovery typically require measuring distances between speech segments [1,2,3]. To avoid the computational expense of the traditional Dynamic Time Warping method, recent papers often use Acoustic word embeddings (AWEs), which represent variablelength segments as fixed-dimensional vectors [4,5]. These can then be compared quickly using measures such as cosine similarity.
An effective AWE algorithm will embed different instances of the same word close together in the vector space, and instances of distinct words further away. One of the main challenges is how to encode the sequential information from the speech signal into a vector space that has no inherent sequential structure. The representation should encode not just which phones are present, but their ordering-so that words like task, stack, cast, and cats, which all contain the same set of phones, will have distinct clusters in the representational space. 1 Various approaches have been developed that use supervision from known word pairs [6,7,8,9], but here we focus on unsupervised learning of AWEs, where only raw audio is available [10,11,12]. This scenario is potentially important for speech applications in lowresource languages.
A common baseline that preserves sequential order while extracting a fixed-dimensional representation is subsampling (see, e.g., [12,13,14]): selecting a fixed number of (usually equally spaced) frames and concatenating them. Subsampling is simple and fast to compute, but, depending on the input frame size and number of samples, it can lead to prohibitively large embeddings and/or loss of phonetic information. Moreover, it does not perform as well on word discrimination tasks as newer learning-based approaches. Unsupervised learning-based methods typically work in two steps: first, apply an unsupervised term detection (UTD) system [15] to identify similar pairs of segments that are likely to be the same word or phrase, then use the pairs as a noisy set of positive examples to train a neural network. Network architectures vary, but the basic idea is to train the system's representations to make the positive examples closer together in the space [11,12] (and in some models, also to separate additional negative example pairs [16,17]).
Though effective, this learning-based approach relies on running UTD on the target language, which itself is computationally intensive and sensitive to differences in input features [12]. Here, we explore whether using newer self-supervised speech representations, available as pre-trained models [18,19], may obviate both the UTD step and the need for specialized models to learn unsupervised AWEs. We hypothesize that the contextualized speech representations learned by these models will implicitly encode the sequential information needed for AWEs (e.g., by capturing within each frame the local acoustic effects of coarticulation, and/or information at a longer timescale that is needed to reconstruct the masked input during pretraining-where average mask span is nearly 300ms [19]). If so, then it should be possible to create effective AWEs with much smaller dimension than subsampling just by using simple pooling methods such as mean-or max-pooling. While pre-trained models and these pooling methods are widely used across many applications, as far as we know this paper is the first to compare and analyze them for creating AWEs.
We evaluate AWEs created using different pooling operations on the representations from two English pre-trained models-HuBERT (HB) and wav2vec 2.0 (W2V2)-and one multilingual pre-trained model (XLSR-53, which also uses the W2V2 architecture). Using a standard word discrimination task, we test on English, Xitsonga, Mandarin, and French-where the latter three better represent a lowresource target language scenario, where a large pre-trained model on that language is unlikely to be available.
In accordance with our hypothesis, we find that on our English test set, AWEs created by mean-pooling the HB representations perform almost as well as the state-of-the-art learned pooling model (MCVAE [20]) with equivalent dimensionality; and outperform subsampling, despite having a much lower dimensionality. Mean-pooled W2V2 representations underperform HB, but are still better than subsampled ones, and considerably better than the MFCC baseline.
Our experiments on other languages show that (1) the multilingual W2V2 representations work better than the monolingual English ones, but still underperform HB (for which only an English model is available); (2) unlike on English, mean-pooling does not outperform subsampling with the HB representations, although it comes close; and (3) when equated on dimensionality, the HB representations are not quite as good as those learned by the best recent models, but perform surprisingly well given that (unlike these models) they require no training on the target language at all.
Overall, our results indicate that the right self-supervised model capture some sequential information needed for AWEs, making simpler pooling methods effective. While these contextualized representations don't generalize fully to other languages, they still work well with no training required.

OVERVIEW OF THE APPROACH
Our approach starts by encoding the corpus using a pre-trained model. The embedding for a given word is created by extracting the model's representations from the start to the end of that word and pooling these to create a fixed dimensional representation.
More formally, let us define two word segments x 1 s 1 :t 1 and x 2 from utterances x 1 and x 2 , where si and ti are the start and end times of the word from x i . We encode both utterances using a contextual self-supervised encoder f , yielding z 1 = f (x 1 ) and z 2 = f (x 2 ). We then pool the encoded representations of each word using a pooling function g, to obtain embeddings c1 = g(z 1 s 1 :t 1 ) and c2 = g(z 2 s 2 :t 2 ). We experiment with four different pooling functions: subsampling, argmax, sum, and mean.
Our encoder models are HB and W2V2, two recent selfsupervised models based on Transformer architectures with latent states. In both models, CNN layers are used to encode the speech signal into audio features, followed by Transformer layers trained using a BERT-like masked language modeling objective (masking some of the input frames). Rather than using the context to predict exactly these frames, both models aim to predict quantized latent units. In W2V2 these are learned jointly in the neural model, while in HB training iterates between a separate clustering step (using K-means) and training the neural network to predict these clusters. Our pooling functions are applied to the frame-level representations from the Transformer layers.
To evaluate the AWEs, we use the same-different word discrimination task (henceforth, same-diff ) [21]. For each pair of embeddings (c1, c2), we measure their cosine similarity and compare this value to a threshold to decide whether both embeddings belong to the same word type. We repeat this process across all possible pairs of a given set of word instances. By varying the similarity threshold across all possible values, we obtain an ROC curve; the final evaluation measure is the Average Precision (AP), or area under this curve.

EXPERIMENTS
We experiment with frame-level representations from three models: English W2V2 and HB, and multilingual W2V2. For the English models, we use Wav2Vec 2.0 Large and HuBERT Large from the official repository 2 , which are both pre-trained on the 60k hour split from the Libri-light dataset [22] and have 317M and 316M parameters respectively. For the multilingual model, we use XLSR-53 from the official repository, which is trained on 53 languages (including Mandarin and French, but not Xitsonga) and has 317M parameters. All models have a contextual representation with 1024 dimensions, so the c's also have 1024 dimensions when g is argmax, sum or mean. For subsampling, we concatenate 10 equally spaced frames, resulting in a 10240-dimensional embedding. Except where noted, frame-level representations, the z's, are normalized by subtracting the mean and dividing by the variance of the evaluated set.
All models have 23 Transformer layers. To save time and computation, we limit our study to layers 1, 11, 15, 19 and 23 and choose the best layer using the development set when available.
We test our AWEs on English, Mandarin, French, and Xitsonga (a low-resource Bantu language spoken in southern Africa). For the English experiments, we focus mainly on a cross-domain setting, testing on the Buckeye corpus of conversational speech [23] (whereas the pre-training corpus, Librispeech, is read speech). We use the dev and test splits (6h each) defined by [24] 3 . We also include some in-domain results from Librispeech [25] dev-clean and test-clean (5.4h each). The Xitsonga data comes from the NCHLT corpus [26], which contains 2.5 hours of read speech and no dev/test split. For Mandarin and French we use the test sets from the ZeroSpeech Challenge 2017 [27] 4 , containing 2.5 and 24 hours of read speech respectively. With the larger French set, we created separate dev and test sets as described below.
To perform the same-diff evaluation, start and end timestamps for each word are needed. The Buckeye and Xitsonga corpora include manually corrected timestamps; for the French and Mandarin corpora we use the Kaldi forced alignments provided; and for Librispeech we obtained alignments using the Montreal Forced Aligner [28]. Following [24], we evaluate using the words from each split that are at least 5 characters and 0.5 seconds long. 5 Since the French data is larger, we created separate dev and test 6 sets from it by randomly sampling 4000 of the relevant word tokens for dev and another 4000 (without replacement) for test. 7 Where possible, we compare to results from [12] and [20], two recent papers on unsupervised acoustic word embeddings who evaluated on English and Xitsonga. [20] appear to have the best published results on these languages using an architecture they call Maximal Sampling Correspondence VAE, while [12] use one of the most well-studied pooling architectures, the CAE-RNN [11], and provide several useful baselines. Both architectures learn from UTD pairs extracted from the target language data (as described in the Introduction), and [12] explore different frame-level features that are used to represent the UTD pairs as input to the CAE-RNN. Their best results are obtained using input features learned using Contrastive Predictive Coding (CPC) [29], a self-supervised model. Thus, those results combine self-supervision with a learned pooling function, whereas we focus on self-supervision alone (but using pre-trained W2V2 and HB instead of training CPC on the target language).

Monolingual setting
We first present results on English, focusing on the HB representations, which we found to work better than W2V2. (Selected W2V2 results are presented for comparison in Section 3.2). Preliminary experiments showed that representations from layer 19 worked best for English, so all results in this section use that layer from that model.  work on related tasks [30,31], we find that normalization greatly improves performance, so we use it in all remaining experiments (including with W2V2, where it also helps). More importantly, we find evidence to support our hypothesis that contextualized representations from self-supervised models implicitly encode sequential information. Specifically, we see that the mean-and sum-pooling strategies work as well as sub-sampling, despite the latter having 10 times more dimensions and explicitly modeling sequential order. Since mean and sum work equally well, we focus on mean-pooling (as compared to subsampling) in the remainder of the paper. So far, our AWEs have at least the same dimensionality as the pretrained model representations, i.e. 1024 dimensions. However, some downstream tasks require fewer dimensions due to computational constraints [13], and many previous AWE systems focus on fewer dimensions. Therefore, for comparison with previous systems, we use PCA to reduce the dimensionality of the frame-level representations to 130 dimensions. For the subsample-pooled AWE, we perform a 13 dimensions PCA to the frame-level representations and then perform subsampling. For mean-pooling, we reduce the frame-level dimensionality to the target dimension (e.g., 130). Figure 1 (bottom) shows that comparison, where we see that for embeddings of the same size down to 130 dimensions, mean pooling always outperforms subsampling, though the benefit is less for smaller embeddings. The figure shows development results, but we also computed test set results with 130 dimensions in order to compare to the best published results for an unsupervised AWE model (MCVAE) [20] (also with 130 dimensions). The AP score of the HB mean-pooled embeddings (35.2%) is close to the MCVAE (39.5%) despite having a much simpler pooling strategy-although we note that the MCVAE is trained using a much smaller amount of English data than the pre-training data for HB.

Cross-lingual setting
Although the results on English are promising, unsupervised AWEs are more likely to be useful for low-resource languages, where a large pre-trained model may not be available for the target language. In this section, we investigate how well AWEs obtained from English or multilingual pre-trained models can work on other languages. We test on two languages that are included in the multilingual pre-training data to differing degrees (Mandarin and French, with 12h and 1686h, respectively, included in the 56k total hours of pre-training data), and and Multiple languages (ML) and mean as pooling mechanism. All results are computed on the development. We use mean pooling, but relative perfromance across layers is consistent with subsampling.
one truly low-resource language (Xitsonga), which is not included in pre-training (nor are any other closely related languages).
We start by determining the best layers to use and whether that differs across languages or models. Figure 2 (top) shows the development set performance on same-diff task across different layers on the English datasets: Librispeech and Buckeye. First, we observe that the result is consistent with previous work on W2V2 [32], which showed that phone and word identities were encoded most strongly in the middle layers. For HB, the trend is different: the later layers generally have better AP scores, suggesting that phone and word identities are better encoded there. In terms of the best performance, HB outperforms both versions of W2V2. Since the English models are pre-trained with the same amount of data and have similar numbers of parameters. This suggests that the difference in performance is due to differences in the pre-training objective.
Next, we repeat the layer-wise analysis on French (dev set) and Xitsonga. Since Xitsonga and Mandarin have no dev sets, we do this analysis on the Xitsonga test set but leave out Mandarin to avoid over-using test sets. As on English, we see that for W2V2 the best performance is near the middle, while for HB the later layers are better-though unlike in English, the final layer appears to be best. We can also see that, for the best-performing layer of each model, multilingual W2V2 works better for these non-English languages than English W2V2, but the improvement is not enough to outperform the English HB model. We hypothesize that a multilingually trained HB model might do even better, but no such model is currently available.
We now compute test set performance on all data sets using the best layers from previous analyses (HB: layer 19 for English, 23 for other languages; W2V2: layer chosen using dev set where available, otherwise layer 11). Figure 3 shows the AP scores for all models, using mean-pooling (1024 dims) and subsampling (10240 dims). For the best HB model, we also include subsampling after PCA to match the dimensions of mean-pooling. Unlike in English, subsampling outperforms mean-pooling for all other languages. This is likely because the English contextual information learned by the model doesn't generalize fully to languages with different phonotactic patterns (see Section 3.3 for analysis); therefore the explicit sequential modeling provided by subsampling can still help. Interestingly, reducing dimensionality prior to subsampling (white bars in Figure 3) hardly reduces the performance on the non-English languages, in contrast to the large drop for English. Therefore mean-pooling and subsampling yield    [12,20], using 130-dimensional embeddings (Table 1). While the crosslingual HB embeddings are slightly worse than previous results using CAE-RNN and MCVAE pooling, they don't require any training on the target language. That said, the CPC features [12] and pooling architectures [12,20] were trained on only a few hours of data, so in future it would be worth comparing CPC and HB representations trained on comparable amounts of data, and/or applying CAE-RNN or MCVAE to HB representations.

Qualitative analysis
We hypothesized that AWEs from self-supervised models would implicitly encode sequential information. Our quantitative results support this hypothesis and suggest that some of this information is language-specific. However, we wish to know whether the results are simply due to good encoding of coarticulation or whether longer-distance sequential information is also encoded. To explore this question, we looked at five words that (according to CMUDict) contain the same set of four phones in different orders. Figure  4 (left) visualizes the mean-pooled HB AWEs for all instances of these words in Librispeech test-clean, dev-clean, and train-clean-100. We observe that the instances are clustered according to their respective word types, providing further evidence that the AWEs encode sequential ordering in some way.
More intriguingly, the most frequent word (asked) forms three sub-clusters, and when we examined them, we found evidence that they encode word-level context rather than just immediate phonetic context (as from coarticulation). To illustrate, Figure 4 (right) color- codes the instances of asked depending on the following word, and shows that instances followed by pronouns (shades of blue) cluster together. Critically, instances followed by me cluster with the other pronouns, rather than with the instances that are followed by mister, even though the latter shares its initial phone with me while the other pronouns do not. Overall, this analysis suggests that the HB AWEs encode much more than just coarticulation, which may explain why mean-pooling is less successful in the cross-lingual setting. It is also unclear whether sub-clustering frequent words by word-level context is desirable for AWEs; this may depend on the task.

CONCLUSION
We hypothesized that self-supervised frame-level representations contain sufficient context so we do not need to model sequential order to construct AWE. Our results on conversational English confirm that this can be true, but depends on the model: for HuBERT (but not W2V2) AWEs constructed using mean-pooling outperform subsample-pooling despite having 10 times fewer dimensions. Hu-BERT representations also perform better overall, suggesting that the nature of the input features has a strong influence on AWE quality, regardless of the pooling method.
When we applied HuBERT to other languages, we found that although the encoded context trained on English does not fully generalize to those languages, mean-pooling worked almost as well as subsampling. With unreduced dimensionality, these methods perform similarly or better than state-of-the-art unsupervised AWE models, though performance degrades when equating dimensionality. Unlike the comparison models, these results are obtained with no training on the target language, so a promising future direction would be to explore whether results on other languages can be improved by adapting the HB model to the target language by continuing self-supervised training on a small amount of unlabeled target language data.