Active Learning for Discovering Complex Phase

Diagrams with Gaussian Processes

Max Zhu1,† , Jian Yao2,† , Marcus Mynatt3 , Hubert Pugzlys3 , Shuyi Li3 , Sergio Bacallado1 ,
Qingyuan Zhao1,* , and Chunjing Jia3,*

1 Statistical
Laboratory, University of Cambridge, Cambridge, CB3 0WA, United Kingdom
arXiv:2409.07042v1 [[Link]-ph] 11 Sep 2024

2 Department of Physics, Southern University of Science and Technology, Shenzhen, 518055, China
3 Department of Physics, University of Florida, 2001 Museum Road, Gainesville FL 32611, United States
† Equal Contribution
* qyzhao@[Link]; chunjing@[Link]

ABSTRACT

We introduce a Bayesian active learning algorithm that efficiently elucidates phase diagrams. Using a novel
acquisition function that assesses both the impact and likelihood of the next observation, the algorithm iteratively
determines the most informative next experiment to conduct and rapidly discerns the phase diagrams with multiple
phases. Comparative studies against existing methods highlight the superior efficiency of our approach. We
demonstrate the algorithm’s practical application through the successful identification of the entire phase diagram
of a spin Hamiltonian with antisymmetric interaction on Honeycomb lattice, using significantly fewer sample
points than traditional grid search methods and a previous method based on support vector machines. Our
algorithm identifies the phase diagram consisting of skyrmion, spiral and polarized phases with error less than 5%
using only 8% of the total possible sample points, in both two-dimensional and three-dimensional phase spaces.
Additionally, our method proves highly efficient in constructing three-dimensional phase diagrams, significantly
reducing computational and experimental costs. Our methodological contributions extend to higher-dimensional
phase diagrams with multiple phases, emphasizing the algorithm’s effectiveness and versatility in handling complex,
multi-phase systems in various dimensions.

1 Introduction and Background

The phase diagram represents an indispensable map for understanding the different states of matter and the phase
transitions between them under varying conditions, such as temperature, pressure, chemical concentration etc. This
conceptual tool is not only fundamental in thermodynamics but also plays a pivotal role in material physics for
the understanding of the interplay of different microscopic degrees of freedom. On the practical side, it aids the
design of functional materials and the exploration of emergent phenomena for next-generation energy materials,
guiding the synthesis and characterization of materials such as alloys [1], perovskites, and strongly correlated
electronic systems [2]. For instance, in the realm of unconventional superconductivity, phase diagrams serve a
critical role in identifying the conditions (temperature, doping concentration, pressure, etc) under which materials
exhibit superconductivity and display the competing phases such as antiferromagnetism, fermi liquid phase, etc.
[3]. In the design of perovskite solar cells, which represent a promising photovoltaic technology, the ability to
design new halide perovskite systems could rely on stabilizing the desired functional phases within specific strain
or pressure ranges. [4, 5] Likewise, in the realm of magnetic materials, comprehending the phase diagram across
various microscopic parameters such as spin exchange and Kitaev interactions, as well as external magnetic field, is
essential for exploring novel magnetic states of matter, including quantum spin liquids, non-collinear magnetism,
and magnetic skyrmions, etc. Thus, the phase diagram serves as an invaluable reference for scientists and engineers
striving to manipulate the physical states of matter toward innovative ends, and efficiently determining the phase
diagram for complex systems across different materials systems is crucial.
Physicists and materials scientists have invested significant effort in developing state-of-the-art experimental
techniques and numerical simulations to obtain the phase for each parameter set. To obtain the entire phase diagram,
the traditional approach, whether experimental or theoretical, typically entails a time-consuming grid search over the
entire phase space. However, this process can become highly inefficient as the dimension of the parameter space and
the number of phases increase, and it can become prohibitively expensive or time-consuming when sophisticated
experiments or simulations are needed to evaluate the phase. To our knowledge, despite much methodological
progress about experimental design in statistics and machine learning, little attention has been given to the problem
of discovering complex phase diagrams until recently [6, 7, 8]. Here we present an active learning framework to
increase sampling efficiency.
Active learning [9] is an iterative process in which an algorithm uses existing data to determine how to sample new
observations. Originally proposed for complex optimization problems [10], active learning has found applications in
a broad range of problems including optimizing material design [11]. Our work is motivated by the active learning
method in Dai and Glotzer [12] that uses a Gaussian Process (GP) regression [13] to prioritize sampling new points
near the predicted phase boundary and of high uncertainty. However, the acquisition function in Dai and Glotzer
[12] is based on a simple heuristic and does not fully use the GP posterior, and their method is limited to modeling
diagrams with just two phases. Another recent work by Tian et al. [14] investigated phase diagrams with more than
two phases using a classifier to predict phases. They propose two methods, one based on Support Vector Machines
that does not quantify the uncertainty of boundaries, and a second method that quantifies uncertainties but requires
direct observation of the phase transition point at the phase boundary, rather than relying on samples taken from
inside each phase. Our proposed method can quantify uncertainties without requiring direct observations of phase
boundaries.
We propose a general Bayesian active learning framework for efficiently discovering phase diagrams. Our
method is particularly adept at addressing phase diagram challenges encountered in physics and materials science,
where each phase is contiguous in parameter space and the occurrence of isolated "islands" of phase within another
phase is improbable. Compared to limited existing approaches, our method has the following features: (1) the
phase diagram can contain any number of phases; (2) the acquisition function is purposefully designed using the GP
posterior to estimate the expected change in phase diagram after sampling a new point; (3) we do not require direct
measurements the phase boundary and allow noisy measurements of the phase at the sampled parameter values. As
a result, our algorithm converges to the true phase diagram considerably faster than existing methods in simulations.
We apply our algorithm to investigate the phase diagram of a spin system with antisymmetric interaction on a
monolayer honeycomb lattice under the influence of a magnetic field. Skyrmions, as an exotic state of matter with
topological spin configuration, hold significant promise for next-generation technological applications in areas such
as data storage, spintronics, sensor technology, and quantum computation, etc.[15, 16] Our algorithm efficiently
identifies the skyrmion phase, which is sandwiched between a spiral phase and a polarized phase, using only a small
sample in the Hamiltonian’s high-dimensional parameter space. Specifically, our algorithm identifies the entire
phase diagram, including skyrmion, spiral, and polarized phases, with an error of less than 5%, using only 8% of the
total possible sample points. The results are robust in both two-dimensional and three-dimensional phase spaces. As
demonstrated by the effective sampling for obtaining the entire phase diagram, we show that our active learning
algorithm has the great potential to revolutionize the time-consuming phase diagram inquiry if it gets embedded
with the automated sampling in phase space for sophisticated experiments or simulations. This integration could
significantly accelerate the discovery of emergent phases and the design of functional materials, particularly those
with high-dimensional tuning parameters or where each experiment or simulation is time-intensive.

2 Methodology
Our objective is to model a phase diagram with n phases (e.g. state of matter) determined by k parameters x ∈ Rk
(e.g. temperature or interaction strength), where the true phase y(x) ∈ {1, ..., n} varies with x. We assume the user
can conduct experiments (or simulations) at a point x to give a possibly noisy measurement ỹ = (p1 , .., pn ) ∈ Rn ,
where pm = p(y = m | x = xt+1 ) is the probability of phase m being the true phase at point x. Our active learning
algorithm takes existing observations T = {xi , ỹi }ti=1 at step t to determine the next point xt+1 to sample, and the
user conducts an experiment at xt+1 to give observation ỹt+1 that is added to the dataset. This process is repeated

2/16
until a phase diagram with satisfactory precision is obtained. A good active learning algorithm would suggest new
points that allow the user to discover the phase diagram with as few measurements as possible.

Observations For noiseless observations, the phase is exactly measured so the measurement posterior consists
of a one-hot vector ỹ(x) = (0, ..., pm = 1, ..., 0) if the measured phase is y = m. For noisy observations, ỹ is the
estimated probability of the true phase, including experimental error. However, real-world experiments often do
not directly return a full distribution. In this case, ỹ is the measurement posterior distribution derived from the
experimental outcome and statistical error. For example, each experiment may return a single phase that is incorrect
with probability ε that can be estimated from prior knowledge, e.g. through analysis of experimental uncertainty or
knowledge about experimental noise levels. If phase m is deemed most probable in this case, the vector ỹ could
consist of pm = 1 − ε and p j = ε/(n − 1), j ̸= m for all other phases. Experimental error can vary with experimental
parameters, ε = ε(x) to allow for increased experimental errors near the phase boundary. This flexible setup allows
the user to incorporate more detailed knowledge of measurement errors, for example, when multiple experiments are
conducted at a single parameter value or if different phases have different measurement errors.

Gaussian Process Regression To predict the true phase at unmeasured parameter values and quantify its
uncertainty, we use a Gaussian Processes (GP) regression model. We give a brief review of GP regression here and a
more detailed introduction can be found in [13].
A Gaussian process (or a Gaussian random field, when the dimension k > 1) is a collection of random variables
{a(x) : x ∈ X } such that for any finite set of indices x1 , . . . , xm ∈ X , the marginal distribution of (a(x1 ), . . . , a(xm ))
is multivariate normal. In GP regression, we observe a set of observed responses a = (a(xx1 ), .., a(xxn )) at points
X = (x1 , . . . , xn ) and we would like to predict the responses a∗ = (a(xn+1 ), . . . , a(xn+n∗ )) at new points X∗ =
(xn+1 , . . . , xn+n∗ ). The GP prior assumes that the responses follow a Gaussian process with covariance kernel
function K, that is, the covariance between two points a(xi ) and a(x j ) is given by Σi j = K(xi , x j ). The prior joint
distribution is (assuming zero prior mean):
     
a 0 Σ Σ∗
∼N , T , (1)
a∗ 0 Σ∗ Σ∗∗

where Σ = K(X, X) is the matrix of the covariances evaluated between all pairs of known points X and Σ∗ = K(X∗ , X)
and Σ∗∗ = K(X∗ , X∗ ) are defined similarly. When a contains noisy measurements of the "true" response, a diagonal
matrix can be added: Σ = K(X, X) + σ 2 I where σ 2 is the noise level. The posterior distribution of a∗ is then given
by

a∗ | a ∼ N (ΣT∗ Σ−1 a , Σ∗∗ − ΣT∗ Σ−1 Σ∗ ). (2)

Since this posterior is normal, the uncertainty of the prediction is quantified by the covariance matrix.
In the application below we choose the Matern 1/2 kernel [17] as the covariance kernel K. This is a stationary
kernel, so the covariance between two points depends only on the distance between them. The kernel has two
parameters—length scale and variance—that are fitted using maximum marginal likelihood estimation (along with
σ 2 if the response is noisy). To avoid cluttering, we will drop the stars in the posterior distributions below.

Our Model We next describe how we fit the GP regression model to phase measurements. We model the probability
distribution of the measured phase diagram by a softmax function of am (x), m ∈ {1...n}, where am (x) is the logit for
phase m at point x. The logits are combined to give the predicted observation probability for each phase m as

exp(am )
p(y = m | x, a1 , . . . , an ) = n , (3)
∑ j=1 exp(a j )

so a higher logit value (relative to other phases) gives a higher probability. The logits a1 (x), . . . , am (x) are assumed
to be independent GPs a priori. This allows us to fit the GP regression for each phase separately as described below.

3/16
Figure 1. Visualization of our model internal state on a phase diagram with four phases. (a) Predicted most likely
phase (shaded) with previous observations (crosses), colored by the phase observed. The recommended point to
sample is shown as an orange dot. (b) Acquisition function rating points to sample. (c) Predicted error (%) of
observing different phases from predicted.

Since our dataset T consists of a series of probabilities for every phase, we transform the probabilities into
logits using the inverse of Equation 3:
n
1
am = log(pm ) − ∑ log(p j ). (4)
n j=1

The logits can be scaled by an arbitrary additive constant, so we choose to normalize them to sum up to 0. This
is consistent with the choice of zero prior mean for the GPs. However, a practical issue with equation (4) is that
noiseless phase measurements with pm ∈ {0, 1} result in infinite logits. To address this, we clip the probabilities
between 0.01/(n − 1) and 0.99.
Next, we fit GP regression models on the calculated logits to give a posterior distribution of the logits using
equation (2). The posterior distribution of the phase at point x is calculated as
Z
p(y = m | x, T ) = p(y = m | x, a1 , ..., an ) ∏ p(a j | x, T )da, (5)
Rn j

so the the uncertainty about the logits is properly aggregated. If all logits have a high variance (i.e., low certainty),
all the probabilities will be around 1/n. The integral in equation 5 is numerically approximated using Gauss-Hermite
quadrature [18] and can be interpreted as a probabilistic prediction of the experimental outcome at point x. The
predicted phase is taken as the phase with the highest predicted probability,
PT (x) = argmax p(y = m | x, T ). (6)
m

Evaluating this over a grid generates a predicted phase diagram, PT . The uncertainty at a given point can be
characterized by the probability that the true phase is anything other than the most likely phase. This predicted phase
diagram PT is taken as the model predictions, while the predicted uncertantity is 1 − argmaxm p(y = m | x, T ), the
probability of measuring any other phase.
Acquisition function In active learning, the next point to sample xt+1 is determined in a greedy fashion by
maximizing the acquisition function A(x) that measures how much additional information a new measurement at x
can provide: xt+1 = argmaxx A(x). A good choice of acquisition function leads to a more efficient experimental
design.
We propose to use expected change in area/volume of the phase diagram as the acquisition function. Consider
the following distance function between two phase diagrams P1 and P2 : (let V be the total volume of the phase space)

1 1
Z
D(P1 , P2 ) = 1(P1 (x) ̸= P2 (x))dx ≈ ∑ 1(P1 (x) ̸= P2 (x)). (7)
V |G | x∈G

4/16
where G is a pre-chosen grid of points used for numerical approximation. Our acquisition function is defined as
A(x) = Eŷ [D(PT , PT ∪{x,ŷ} )], (8)
the expected distance between the current most probable phase diagram PT and the most probable phase diagram
after a "fantasy experiment" is performed at point x. The expectation in Equation (8) is taken over the posterior
distribution of the "fantasy measurement" ŷ at x given by the GP regression model. "Fantasy experiments" are
hypothetical observations that are distributed according to our current model (Equation (5)) that simulate the outcome
of real experiments,

ŷ = (ε̂/n, ..., pm = 1 − ε̂, ..., ε̂/n), m ∼ p(y|x, T ), (9)

where ε̂ is a user chosen hyperparameter representing a rough estimate of experimental noise and the main measured
phase j is distributed according to our model. Ideally, ε̂ would equal the real experimental noise, but this is unknown
so we set it to a constant using prior knowledge about experimental noise. For the noiseless case, we choose ε̂ = 0.01
because setting it too low may result in numerical instability. Fantasy observations mirror the setup described in the
Observations section.
To compute A(x), we fit our model on our current observations T and generate the predicted phase diagram
PT using Equation 6. Then, for each phase m, we fit a new model on the current dataset with an extra fantasy
observation, {x, ŷ = ŷm }, where ŷm is the fantasy observation probabilities for phase m (Equation 9), and generate
predictions giving PT ∪{x,ŷm } . Finally, A(x) is computed by taking the mean distance between the current phase
diagram and each possible fantasy phase diagram, weighted by the probability that each fantasy diagram occurs.
This acquisition function finds points that will maximally change the phase diagram, by combining the impact of
an observation on the entire phase diagram with the probability of that observation and balances the exploration
of unknown regions with the refinement of known boundaries. Compared to other commonly used acquisition
functions, such as KL divergence or Expected Error Reduction [19], our acquisition function prioritizes discovering
the phase boundary while placing less emphasis on points far away from the phase boundary. The method in Dai and
Glotzer [12] differs from ours as it finds the points that are locally maximally uncertain, instead of considering the
impact of an observation on a region surrounding the point. This process is visualized in Figure 1. Our acquisition
function is shown in Figure 1(b), while Dai and Glotzer would pick a point that maximizes the uncertainty in Figure
1(c).
Implementation Details We summarize the processes of our proposed active learning approach for phase diagram
identification in algorithmic flowchart as shown in Figure 2. We implement our method using PyTorch [20].
Appendix 7.2 describes implementation details and several numerical optimizations to speed up the algorithm.
With these optimizations, each step of active learning requires about 25 seconds on an 8-core processor for a
three-dimensional phase diagram.

3 Synthetic Experiments
Our model is evaluated against Dai and Glotzer [12] (Gaussian process regression or GPR)
and a modified version of Approach 1 in Tian et al. [14] (support vector machine classification or SVC)
as baselines on two synthetic examples. Like our method, GPR is a Bayesian Gaussian Process method but is
only suitable for binary classification and cannot directly handle experimental uncertainty. The SVC method uses a
support vector machine to predict classification boundaries and places observations farthest from other observations
on the boundary. The SVC cannot account for uncertainty, quantify prediction uncertainty or discover unknown
isolated phases since the model will not search regions not on the boundary. Further details of these baselines are
given in Appendix 7.1. These experiments are restricted to k = 2 for easier analysis and visualization. The first
experiment is on a noisy phase diagram, while the second has several phases.
Two phase problem: Firstly, we evaluated the models on a noisy toy phase diagram of two phases, P(x1 , x2 ) =
1(x2 > sin(0.5πx1 )) where x1 , x2 ∈ [−2, 2] represent experimental parameters. Uniform noise, ε ∼ [−0.2, 0.2] was

5/16
 1. Generate phase diagram from 2. Calculate acquisition 3. Select best point for 4. Update predictions, and
current observation dataset function for every point experimentation and conduct repeat process
experiment

Select one point as

an example Fantasy Observations to compute

Simulate Phase 1 2 3

Outcome

Probability 0.48 0.22 0.30

Distance 0.029 0.057 0.050

Acqusition 0.042

Figure 2. A schematic flowchart of the algorithm for our proposed active learning approach to phase diagram
problems.

Figure 3. (a) A phase diagram of four phases used for testing, with phases shown with shaded color. Predicted
phase diagram after 100 observations from our model (b) and SVC (c), with predicted phase shaded and
observations as crosses colored by observed phase. Note the SVC repeats a large number of observations and does
not attempt to search the green region for undiscovered phases.

6/16
 a b

Figure 4. Fractional error between predicted and true phase diagrams versus number of samples taken for phase
diagrams of (a) two phases and (b) four phases. Shaded area in (a) shows standard deviation over 6 runs with
different random noise, while (b) is noiseless.

added to x1 , x2 every observation to simulate experimental errors. Models were initialized with 2 observations. We
assume measurement error is known for this experiment, so we input the error probability into our model, while the
baselines cannot use the measurement error. Figure 4 and Table 1 show the percent difference between the predicted
and true phases over the entire phase space (fractional error) vs. the number of observations, averaged over 6 runs.
Our method surpasses the baselines in accuracy after 12 observations and achieves 60% lower final error.
Four phase problem: Our model and Tian (SVC) were evaluated on a more complex phase diagram of four
spaces inspired by the pressure-temperature phase diagram of materials (Figure 3 (a)). Measurements are noiseless.
Both models were initialized with 2 observations, so there are 2 phases hidden at the start so the model must discover
the location of the two hidden phases. Figure 4 (b) shows the error vs. number of observations and final predictions
are shown in Figure 3 (b) and (c). Our model substantially outperforms the SVC after 22 steps, and perfectly predicts
the phase diagram after 83 steps. SVC predictions are unaligned with observations, which was also observed in
Ref [14]. Figure 1 visualizes our model’s predicted phase diagram, acquisition function, and error probability. Figure
A1 in the Appendix shows the evolution of our model’s predictions.

4 Application on Phase Diagram of Magnetic Materials: Heisenberg Hamiltonian with

Antisymmetric Interaction
A motivating application is the exploration of exotic phases in frustrated magnetic systems. One particular phase
that has recently garnered significant attention is magnetic skyrmions. In the two-dimensional scenario, magnetic
skyrmions manifest as real-space topological structures, formed by localized spin vectors that smoothly transition
from their cores to their peripheries while maintaining fixed helicities. The real-space topology of a skyrmion is
1 R
characterized by its skyrmion number Q = 4π m · ( ∂∂mx × ∂∂my )dxdy, where m is the unit vector of local magnetization
and the integral is taken over a two-dimensional space.[21] To investigate the skyrmion phase, we focus on a
Heisenberg model with antisymmetric interaction on the honeycomb lattice:
H = −J ∑ S i · S j − ∑ d i, j · (SS i × S j ) − D ∑(Siz )2 − h ∑ Siz , (10)
⟨i, j⟩ ⟨i, j⟩ i i

where S i represents the classical spin vector at site i with |SS i | = 1 and Siz is z-component. J > 0 is the nearest
neighbor (NN) ferromagnetic coupling, and d i, j = d (zz × u i, j ) is the NN Dzyaloshinskii–Moriya (DM) interaction
[22, 23], where z and u i, j are unit vectors, respectively perpendicular to the magnetic layer plane and pointing from
site i to site j, as shown in Fig. 5 (a). D > 0 represents easy-axis single-ion anisotropy, and h > 0 is the magnetic
field. Besides the skyrmion phase, we have identified another two distinct magnetic phases in this system: the spiral
and polarized phases, in which the spin configuration has magnetic wave vector q and 0 . An example of the spiral
phase and skyrmion phase is illustrated in Fig. 5 (b) (c).

7/16
a b c

d e

h/J

D/J

Figure 5. Phase diagram of Heisenberg model with antisymmetric interaction using our active learning approach.
(a) A schematic of the model Hamiltonian and the microscopic spin configuration for skyrmion. J and DM represent
the spin exchange and DM vector for Dzyaloshinskii–Moriya interaction respectively. Single-ion anisotropy D and
magnetic field h are not shown. (b) and (c) are real-space spin configurations of (b) spiral and (c) skyrmion phases.
Color denotes Sz component (colorbar shown) and arrows denote the spin direction in (a) and (Sx , Sy ) components in
(b) and (c). Here (b) and (c) are obtained at (D, h)/J = (0.056, 0) and (0.028, 0.7), respectively. (d) Phase diagrams
as functions of D/J and h/J, for d = 0.5J and temperature T = 10−4 J, obtained after t observations. Crosses with
colors represent points sampled: blue for the spiral phase, white for the skyrmion phase, and red for polarized phase.
We use purple, green, and yellow to fill out inferred regions of the three different phases respectively. As a
comparison, black lines show the ground-truth phase boundary. (e) Fractional error as a function of the number of
observations (starting from t = 6, with the same initial points) for SVC method (blue), grid sampling (orange), and
our method (green).

In order to determine the magnetic state configuration, we perform atomistic spin dynamics simulation by
VAMPIRE [24] under fixed J = 1 meV, d/J = 0.5 and temperature lowered to T /J = 10−4 . The system size is
chosen to be Lx = Ly = 50a with periodic boundary conditions, where a is the lattice constant. In constructing the
phase diagram for D/J ∈ [0, 0.5] and h/J ∈ [0, 0.25] with 21 × 21 discrete possible points using our algorithm, we
start with 6 random initial starting points, then increase the number of observations to t = 16, 26, 36, 46, 56, 66, 76.
The resulting phase diagrams are illustrated in Fig. 5 (d). As a comparison, the black lines in each plot represent the
true phase boundaries obtained by grid sampling. The fractional error shown in Fig. 5 (e) is calculated by comparison
with the true diagram via 21 × 21 grid sampling. The fractional error decreases as the number of observations grows,
becoming smaller than 5% after t = 37 and smaller than 2.5% after t = 60, as illustrated in Fig. 5 (e). After t = 141,
the fractional error becomes zero, suggesting that our active learning approach now gives a perfect prediction for the
entire phase diagram on the 21 × 21 grid in the two-dimensional phase space. Compared with direct grid sampling
and the SVC method, our algorithm is much more efficient. Grid sampling is implemented by sampling points in
order (sequentially sampling each parameter in order), starting from a half resolution grid before increasing to the
full-resolution grid. Details are given in Appendix 7.1.

8/16
 a b
spiral
skyrmion
polarized

h/J

D/J

d(DMI)/J

c d e

h/J

D/J D/J D/J

d(DMI)/J d(DMI)/J d(DMI)/J

Figure 6. 3D phase diagram of Heisenberg model with antisymmetric interaction using our active learning
approach. (a) The regions of red, green, and cyan colors correspond to the spiral, skyrmion, and polarized phases,
respectively. We also mark sampled observations with magenta, yellow, and blue little circles to represent the three
phases, respectively. (b) Fractional error as a function of the number of observations (starting from t = 6, with the
same initial points) for SVC method (blue), grid sampling (orange), our method (green). (c) (d) (e) 3D phase
diagrams that are highlighted at different values of D.

9/16
 When more than two parameters are considered, the number of possible observations becomes enormous, making
it very difficult to create a phase diagram using grid sampling. To further demonstrate the efficiency of our algorithm,
we use it to generate a three-dimensional phase diagram for the Heisenberg model with varying antisymmetric
interaction d as mentioned above. Instead of setting d/J = 0.5, we now vary d/J in the range of [0, 1], D/J within
[0, 0.5], and h/J within [0.025, 0.25] with 11 × 11 × 10 discrete possible points.
When h = 0, the Hamiltonian exhibits at least two-fold degenerate ground states due to time-reversal symmetry.
The classical spin simulation method we use can result in configurations with multiple domains, complicating the
phase determination for polarized and spiral phases. Therefore, we exclude the case of h = 0 in the 3D phase
diagram analysis. The resulting phase diagrams are illustrated in Fig. 6 (a). The fractional error shown by Fig. 6
(b) green line is calculated by comparison with ground truth phase diagram via 11 × 11 × 10 grid sampling. After
around t = 99, the fractional error becomes smaller than 5%, it converges to zero at around t = 440. The ratio of
number of observations and number of total possible sample points is around 8% to reach the fractional error smaller
than 5% for both 2D and 3D sampling. Compared to direct grid sampling and the SVC method, our algorithm
demonstrates superior performance in both the rate of error reduction and the final error achieved.
To investigate the change of magnetic phases with varying antisymmetric interaction, single-ion anisotropy and
external field, we snapshot a few 2D cuts of the 3D phase diagram in Fig. 6 (a) at easy-axis single-ion anisotropy
D/J = 0, 0.25, 0.5, as shown in Fig. 6 (c) (d) (e). Let us first focus on single-ion anisotropy D = 0, a larger magnetic
field h favors the polarized phase, whereas stronger antisymmetric interaction d causes the spins to form spiral
structures. In the intermediate region, the competition between the magnetic field h and antisymmetric interactions
d could stabilize the skyrmion phase. As the single-ion anisotropy D increases, the phase boundaries shift to the
right, corresponding to higher antisymmetric interaction d, and the polarized phase becomes more competitive.
This means that the polarized phase can be stabilized at higher values of d with increasing D. This behavior is
expected because an easy-axis D favors the alignment of spins in parallel directions (or antiparallel directions for
other systems), requiring a larger antisymmetric interaction d to disrupt the polarized phase. Concurrently, a larger
antisymmetric interaction d is needed to stabilize skyrmion phases as the single-ion anisotropy D increases.

5 Conclusion and Discussion

Our proposed method is able to efficiently determine phase diagrams by identifying where to sample new observa-
tions, using a Bayesian approach to quantify the uncertainty. Demonstrated to be effective for synthetic experiments
and for skyrmion phase diagram problem in both two-dimensional and three-dimensional phase space, this approach
holds promise to significantly accelerate the identification of phase diagrams across a wide array of materials
physics problems, even as the dimensionality of the parameters space increases. Our algorithm is shown to converge
faster than grid sampling and other existing active learning methods. Additionally, our method can adapt to noisy
measurements, making it particularly useful for experimental measurements where empirical errors are common.
We provide our method as a package with source code and a graphical user interface for ease of use.
For future work, we would like to extend our approach to encompass more detailed prior information about
the phase diagram. Users may have prior information about the shape, connectivity or other information about the
phase diagram from before experiemnts are conducted, such as from prior experiments or theoretical work. In these
cases, our Gaussain Process assumption that, a priori, phases are isotropically and uniformly Gaussian distributed
would no longer hold true. Incorporating this information should further increase the efficiency of the active learning
algorithm.
To summarize, our active learning algorithm, if integrated with experimental or numerical methods for the
automated sampling of phase diagrams, facilitates a higher level of automation with strategic and effective sampling
in the phase space. We expect that this advancement will revolutionize scientists’ and engineers’ search for new
materials and states of matter in high dimensional parameters space, allowing for more strategic exploration and
unprecedented efficiency.

10/16
6 Acknowledgments
We thank Sathya Chitturi for the insightful discussion to the project. M. Zhu and Q. Zhao acknowledge support
for this work from GSK, the Cambridge Centre for AI and Medicine, and the Engineering and Physical Sciences
Research Council (Grant EP/V049968/1). C. J. Jia acknowledges the support from Department of Energy, Office of
Science, Basic Energy Sciences, Materials Sciences and Engineering Division, under Contract No. DE-SC0022216,
as well as the support from University of Florida Seed Grant. M. Mynatt acknowledges the support from the
University Scholars Program at University of Florida. H. Pugzlys acknowledges the support from Christopher B.
Schaffer undergraduate research scholarship.

References
1. H. Okamoto, M.E. Schlesinger, and E.M. Mueller. Alloy Phase Diagrams. ASM International, 04 2016. ISBN
978-1-62708-163-4. doi: 10.31399/[Link].v03.9781627081634. URL [Link]
[Link].v03.9781627081634.
2. Elbio Dagotto. Correlated electrons in high-temperature superconductors. Rev. Mod. Phys., 66:763–840, Jul 1994.
doi: 10.1103/RevModPhys.66.763. URL [Link]
763.
3. Callum J. Court and Jacqueline M. Cole. Magnetic and superconducting phase diagrams and transition
temperatures predicted using text mining and machine learning. npj Computational Materials, 6(1):18,
Mar 2020. ISSN 2057-3960. doi: 10.1038/s41524-020-0287-8. URL [Link]
s41524-020-0287-8.
4. Julian A. Steele, Handong Jin, Iurii Dovgaliuk, Robert F. Berger, Tom Braeckevelt, Haifeng Yuan, Cristina
Martin, Eduardo Solano, Kurt Lejaeghere, Sven M. J. Rogge, Charlotte Notebaert, Wouter Vandezande, Kris P. F.
Janssen, Bart Goderis, Elke Debroye, Ya-Kun Wang, Yitong Dong, Dongxin Ma, Makhsud Saidaminov, Hairen
Tan, Zhenghong Lu, Vadim Dyadkin, Dmitry Chernyshov, Veronique Van Speybroeck, Edward H. Sargent,
Johan Hofkens, and Maarten B. J. Roeffaers. Thermal unequilibrium of strained black cspbi₃ thin
films. Science, 365(6454):679–684, 2019. doi: 10.1126/science.aax3878. URL [Link]
org/doi/abs/10.1126/science.aax3878.
5. Feng Ke, Chenxu Wang, Chunjing Jia, Nathan R. Wolf, Jiejuan Yan, Shanyuan Niu, Thomas P. Devereaux,
Hemamala I. Karunadasa, Wendy L. Mao, and Yu Lin. Preserving a robust cspbi3 perovskite phase via pressure-
directed octahedral tilt. Nature Communications, 12(461):461, 2021. doi: 10.1038/s41467-020-20745-5. URL
[Link]
6. Guillaume Deffrennes, Kei Terayama, Taichi Abe, and Ryo Tamura. A machine learning–based classification
approach for phase diagram prediction. Materials & Design, 215:110497, 2022. ISSN 0264-1275. doi:
[Link] URL [Link]
article/pii/S0264127522001186.
7. J. Lund, H. Wang, R. D. Braatz, and R. E. García. Machine learning of phase diagrams. Mater. Adv., 3:
8485–8497, 2022. doi: 10.1039/D2MA00524G. URL [Link]
8. M. Link, K. Gao, A. Kell, M. Breyer, D. Eberz, B. Rauf, and M. Köhl. Machine learning the phase diagram of a
strongly interacting fermi gas. Phys. Rev. Lett., 130:203401, May 2023. doi: 10.1103/PhysRevLett.130.203401.
URL [Link]
9. Burr Settles. Active learning literature survey. University of Wisconsin-Madison Department of Computer
Sciences, 2009.
10. H. J. Kushner. A new method of locating the maximum point of an arbitrary multipeak curve in the presence of
noise. J. Basic Eng., 1964.
11. Turab Lookman, Prasanna V. Balachandran, Dezhen Xue, and Ruihao Yuan. Active learning in materials science
with emphasis on adaptive sampling using uncertainties for targeted design. npj Computational Materials, 2019.
12. Chengyu Dai and Sharon C. Glotzer. Efficient phase diagram sampling by active learning. The Journal of

11/16
 Physical Chemistry B, 124(7):1275–1284, 2020. doi: 10.1021/[Link].9b09202. URL [Link]
10.1021/[Link].9b09202. PMID: 31964140.
13. Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian processes for machine learning. The MIT
Press, 2006. ISBN 0-262-18253-X, 2006.
14. Yuan Tian, Ruihao Yuan, Dezhen Xue, Yumei Zhou, Yunfan Wang, Xiangdong Ding, Liu Juan, and Turab
Lookman. Determining multi-component phase diagrams with desired characteristics using active learning.
Advanced Science, 8:2003165, 11 2020. doi: 10.1002/advs.202003165.
15. Giovanni Finocchio, Felix Büttner, Riccardo Tomasello, Mario Carpentieri, and Mathias Kläui. Magnetic
skyrmions: from fundamental to applications. Journal of Physics D: Applied Physics, 49(42):423001, sep 2016.
doi: 10.1088/0022-3727/49/42/423001. URL [Link]
423001.
16. Kang Wang, Vineetha Bheemarasetty, Junhang Duan, Shiyu Zhou, and Gang Xiao. Fundamental physics and
applications of skyrmions: A review. Journal of Magnetism and Magnetic Materials, 563:169905, 2022. ISSN
0304-8853. doi: [Link] URL [Link]
com/science/article/pii/S0304885322007909.
17. William Kleiber Tilmann Gneiting and Martin Schlather. Matérn cross-covariance functions for multivariate
random fields. Journal of the American Statistical Association, 105(491):1167–1177, 2010. doi: 10.1198/jasa.
2010.tm09420.
18. Peter Jackel. A note on multivariate gauss-hermite quadrature. J. Basic Eng., 2005.
19. Nicholas Roy and Andrew McCallum. Toward optimal active learning through sampling estimation
of error reduction. In International Conference on Machine Learning, 2001. URL [Link]
[Link]/CorpusID:14949756.
20. Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen,
Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito,
Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.
Pytorch: An imperative style, high-performance deep learning library. In Advances in Neural Information
Processing Systems 32. 2019.
21. Stefan Heinze, Kirsten von Bergmann, Matthias Menzel, Jens Brede, André Kubetzka, Roland Wiesendanger,
Gustav Bihlmayer, and Stefan Blügel. Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions.
Nature Physics, 7:713, 2011. doi: 10.1038/nphys2045. URL [Link]
22. I. Dzyaloshinsky. A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics. Journal
of Physics and Chemistry of Solids, 4(4):241–255, 1958. ISSN 0022-3697. doi: [Link]
1016/0022-3697(58)90076-3. URL [Link]
0022369758900763.
23. Tôru Moriya. Anisotropic superexchange interaction and weak ferromagnetism. Phys. Rev., 120:91–98, Oct 1960.
doi: 10.1103/PhysRev.120.91. URL [Link]
24. R F L Evans, W J Fan, P Chureemart, T A Ostler, M O A Ellis, and R W Chantrell. Atomistic spin model
simulations of magnetic nanomaterials. Journal of Physics: Condensed Matter, 26(10):103202, feb 2014.
doi: 10.1088/0953-8984/26/10/103202. URL [Link]
103202.
25. B. Berg and M. Lüscher. Definition and statistical distributions of a topological number in the lat-
tice o(3) σ -model. Nuclear Physics B, 190(2):412–424, 1981. ISSN 0550-3213. doi: [Link]
10.1016/0550-3213(81)90568-X. URL [Link]
pii/055032138190568X.
26. Gen Yin, Yufan Li, Lingyao Kong, Roger K. Lake, C. L. Chien, and Jiadong Zang. Topological charge analysis
of ultrafast single skyrmion creation. Phys. Rev. B, 93:174403, May 2016. doi: 10.1103/PhysRevB.93.174403.
URL [Link]
27. M. Lakshmanan. The fascinating world of the landau–lifshitz–gilbert equation: an overview. Philosophical
Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 369(1939):1280–1300,

12/16
 2011. doi: 10.1098/rsta.2010.0319. URL [Link]
1098/rsta.2010.0319.

7 APPENDIX

Number of observations
Accuracy(%) 10 20 30 40 50
Ours 0.071 0.038 0.028 0.021 0.018
GPR 0.139 0.055 0.051 0.052 0.033
SVC 0.071 0.055 0.043 0.047 0.052

Table 1. Accuracy of predicted phase diagrams for binary classification after each algorithm after different numbers
of observations. Results are averaged over 6 runs. The lowest error for each step is highlighted in bold.

7.1 Details of Baseline Algorithms

In this section, we provide a detailed description of the three baseline algorithms. Dai and Glotzer [12] fit a Gaussian
process p(GP) is fit to binary observations y = {−1, +1} and a heuristic acquisition function rates points, x, to select,
A(x) = var(y)/(|E(y)| + ε) where var(y), E(y) and ε are the GP’s variance, mean and a numerical stability factor.
This acquisition function balances selecting points close to the currently predicted phase boundary where E(y) ≈ 0
and regions of high uncertainty where var(y) is large. The numerical stability factor prevents dividing by 0. Similar
to our method, at every step, A(x) is evaluated and the point with the highest value is selected as the next observation.
Tian et al. [14] presents two methods. Method 2 is able to learn phase diagrams of more than two phases,
however, this is not applicable to our setup since we cannot measure the location of phase transitions directly. Instead,
we use a modified version of Method 1. A support vector classifier is fit to the observations of position and phase.
Phase boundaries are predicted from classification boundaries. For active learning, we sample new points along the
phase boundary that have the highest exponentially weighted average distance from existing samples, which acts as
a proxy for regions of the highest uncertainty. This is because support vector classifiers cannot directly quantify
prediction uncertainty, so we use distance from other observations as a proxy for uncertainty. We hand-tuned the
exponential weighting factor for the best performance. Furthermore, since this model only places observations on
the predicted boundary, it cannot discover disconnected phases that have not been observed.
Grid sampling was implemented by sequentially sampling points starting from the origin, one line of points at a
time. For example, in 2D, the column at x = 0 is sampled, followed by x = 1, for each x coordinate. In 3D, the z
column x = 0, y = 0 is sampled, followed by x = 0, y = 1, for every y column, after which the column x = 1, y = 0 is
sampled until every point has been sampled. Since every point is sampled, this method generates a perfect phase
diagram in the noiseless case. A downside of this grid sampling method is slow convergence since sampling starts
from one direction. To mitigate this, instead of directly sampling from the full-resolution grid, the grid sampling
algorithm first samples from a half-resolution grid and moves to the full-resolution grid after the half-resolution grid
sampling is finished. In the 2D case (Figure 5), the full grid is 21 × 21 so the starting grid was 11 × 11 and for the
3D case (Figure 6), the full grid is 11 × 11 × 10 so the starting grid was 6 × 6 × 4. In both cases, the slowdown when
switching from the low-resolution grid to the full-resolution grid is clear. Nearest-neighbour interpolation is used to
make predictions for incomplete grids.

7.2 Implementation Details

Model details: All of our Gaussian processes, p(ai |x), use the stationary Matern 1/2 kernel, K(x, x′ ) = K(r) =
κ 2 exp(r/ρ), where r is the Euclidian distance between x and x’. The kernel parameters are κ, the magnitude of
correlation between observations, and ρ, setting the lengthscale two observations are correlated. These are fitted

13/16
Figure A1. Observations and predictions of our model on the synthetic phase diagram of four phases after different
numbers of observations. Observations are shown by crosses. The color indicates the observed phase.

using maximum likelihood estimation. We force the fitted kernel and noise parameters to remain positive and also
clip them from becoming excessively large or small to stop numerical instabilities. The Matern 1/2 kernel generates
predictions are continuous, in contrast with the Matern 3/2 or 5/2 kernels which also have (mean square) bounded
derivatives. We chose the Matern 1/2 kernel since phase diagrams are discontinuous and do not vary smoothly. For
further details about Gaussian Processes kernels, we recommend Rasmussen and Williams [13] Chapter 4.
Our acquisition function assumes fantasy observations of phase that are distributed according to current model
predictions (Equation 8). We require ε̂, an estimate of the experimental error for fantasy observations. The true
experimental error is unknown, so users can make a reasonable guess for the experimental error, e.g. based on
previous observations. For our experiments, this is set to a constant ε = 0.01 for noiseless phase diagrams (4 phase
diagram in Section 3 and Skyrmion phase diagram) or ε = 0.10 for noisy phase diagrams. Furthermore, for noiseless
experiments, our algorithm is restricted to not search seen points again while noisy observations can be resampled
multiple times.
Numerical Optimisations: To speed up our method, we make several numerical optimizations. Firstly, sampling
probabilities p(y|x, T ) (Equation 5) requires integrating Equation 3 with respect to n Gaussians, for each of the n
phases giving n n-D integrals. Since probabilities sum to 1, this can be reduced to n − 1 integrals. We efficiently and
deterministically compute the integrals using the Gauss-Hermite quadrature instead of random sampling. This step
is the most time-consuming, so it is performed in parallel using multi-processing.
To find the next point xt+1 = argmaxx A(x), we sample points on a regular grid. Optimizations can be made to
speed up the search. After the current probabilities of observing each phase are computed, we only need to sample
A(x) at points where there is significant uncertainty since very low probability observations are unlikely to impact
the search. We set A(x) = 0 if the probability of observing any single phase is greater than 90%. Also, since we
are using a local kernel for our GPs, the effect of new observations is local so the distance D(PT , PT ∪{x,ŷ(x)} ) is
only sampled for local points within a set distance of x instead of the full area. This is set to be a grid of points with
Euclidean distance within half the search range.

7.3 The criterion to determine phases of the honeycomb lattice

The phases under different parameter sets can be easily identified by observing real-space spin textures (fig). We
have also calculated some physical quantities to precisely identify the different phases. We employ the discrete
version of topological charge Q that identifies the existence of skyrmion phase. This version was first introduced
in [25] and used in [26]. The calculation of Q starts by triangulating the entire lattice and then counting the solid

14/16
 Phase Q M
Spiral low or 0 low
Skyrmion high high
Polarized 0 1

Table 2. Criterion used to assist the determination of three phases

angels Ω∆ for each triangle ∆(SS 1 , S 2 , S 3 ) determined by

Ω∆
exp(i ) = ρ −1 [1 + S 1 · S 2 + S 2 · S 3 + S 3 · S 1 + iSS 1 · (SS 2 × S 3 )], (11)
2

where ρ = [2(1 + S 1 · S 2 )(1 + S 2 · S 3 )(1 + S 3 · S 1 )]1/2 , then Q is given by summing over all of the triangles:

1
Q= ∑ Ω∆ . (12)
4π ∆

Here the solid angle Ω∆ ranges from −2π to 2π and the branch cut is the negative real axis when calculating Ω∆ .
We also calculate standard magnetization:
1
M= Siz , (13)
N∑i

where N is the number of sites in the entire lattice.

In Table. 2 , we show the criterion we use to assist the determination of phases.

7.4 Phase diagram obtained by grid sampling

In Fig. A2, we show phase diagram of the honeycomb lattice obtained by 21 × 21 grid sampling and the phase
boundaries (smoothed). To simulate the spin configuration of the honeycomb lattice, we simulate the atomistic spin
dynamics using in Vampire [24]. Starting from a relatively high temperature, we simulate the evolution of spins
using the Landau–Lifshitz–Gilbert (LLG) equation [27], cooling to set temperature we set with custom timesteps. In
this work, we cool the system from T = 20K to T = 10−4 J(0.00116K) (Fig. A3), at each temperature during cooling
process, the system evolves 104 steps with each step 10−16 s, total 4 × 106 steps in the whole process of a simulation.
For the 3D case, we do 11 × 11 × 11 grid sampling, and we cool the system from T = 30K to T = 10−4 J(0.00116K)
when sampling.

8 Additional information
8.1 Data Availability
We provide code for our method at [Link] This contains
our algorithm and a graphical user interface written in Python for users, as well as code required to reproduce the
synthetic experiments in Section 3.

8.2 Competing Interests

The authors declare that they have no competing interests.

15/16
Figure A2. Phase diagrams of the honeycomb lattice as a functon of single-ion anisotropy D/J (x-axis) and the
magnetic field h/J (y-axis) for fixed Dzyaloshinskii–Moriya interaction d/J = 1/2 and temperature T = 10−4 J
obtained by grid sampling. Black: spiral phase; red: skyrmion phase; pink: polarized phase. Green lines represent
the smoothed phase boundaries.

𝑇𝑇/𝐾𝐾

𝑡𝑡/𝑠𝑠

Figure A3. Temperature as a function of simulation time in our atomistic spin dynamics simulation.

16/16