Unsupervised Learning of Rydberg Atom Array Phase Diagram with Siamese Neural Networks

The Ising model on the square lattice is a simple, well studied, and exactly solvable model from statistical physics that exhibits a phase transition. Thus, it provides the ideal starting point for benchmarking the performance of a phase recognition algorithm. Its Hamiltonian is

$\displaystyle H(S)=-J\sum_{<ij>_{nn}}s_{i}s_{j}+h\sum_{i}s_{i},$

(1)

a function of spin configurations $S=(s_{1},...,s_{n})$. In the following, we focus on the ferromagnetic Ising model $J=1$ and set the external field $h=0$.

SNN phase recognition is an algorithm that is able to detect multiple phases in an unsupervised manner. Hence, we trivially combine two Ising models by overlaying them on the same lattice where we can tune each temperature $T,\tilde{T}$, or in other words the ratios $J/k_{B}T$ and $\tilde{J}/k_{B}\tilde{T}$, independently. The combined Hamiltonian

$\displaystyle H(S)=-J\sum_{<ij>_{nn}}s_{i}s_{j}-\tilde{J}\sum_{<ij>_{nn}}\tilde{s}_{i}\tilde{s}_{j}$

(2)

acts on lattices containing two spins per site, $S=((s_{1},\tilde{s_{1}}),...,(s_{n},\tilde{s_{n}}))$. Since there is no interaction between them, phase transitions trivially occur at lines $(T,\tilde{T})\approx(\ \cdot\ ,2.269)$ and $(T,\tilde{T})\approx(2.269,\ \cdot\ )$ in the phase diagram.

A common Hamiltonian that can be implemented by Rydberg arrays has a form similar to that of a Transverse Field Ising Model, meaning it is sign-problem free and thus amenable to simulation on classical computers. The Hamiltonian acts on a collection of atoms which individually act like 2-level systems, having a ground-state $\ket{g}\equiv\ket{0}$ and an excited, so-called Rydberg state $\ket{r}\equiv\ket{1}$. The atoms are subject to a long-range interaction which is described by a van der Waals (vdW) interaction of the form $V_{ij}\sim\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}}^{-6}$ that penalizes atoms that are simultaneously in the Rydberg state [74]. Additionally, the atoms are subject to coherent laser fields: a detuning $\delta$ which acts like a chemical potential, driving atoms into their Rydberg states, and a Rabi oscillation with frequency $\Omega$ which excites ground-state atoms and de-excites atoms in Rydberg states.

$$H=\sum_{i<j}V_{ij}n_{i}n_{j}-\delta\sum_{i}n_{i}+\frac{\Omega}{2}\sum_{i}\sigma_{i}^{x}$$

(3)

where $n_{i}=\outerproduct{1}{1}_{i}$ is the occupation/number operator, and $\sigma_{i}^{x}=\outerproduct{1}{0}_{i}+\outerproduct{0}{1}_{i}$. The interaction strength is typically parametrized in terms of a Rydberg blockade radius $R_{b}$, which describes an effective radius within which two simultaneous Rydberg excitations are heavily penalized: $V_{ij}=\Omega(R_{b}/\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})^{6}$ [75, 76].

The well-known single-spin-flip Metropolis algorithm is used to generate importance sampled Monte Carlo configurations of the Ising model on a square lattice of size $20\times 20$ with periodic boundary conditions [77]. After initializing a random lattice, we evolve the simulation for 7168 MC steps between drawing samples. It is important to note that neural networks can pick up on any residual correlations, thus relying on conventional auto-correlation measures to determine the independence of lattice configurations is not enough. We produce 92 independent configurations at each of 100 temperatures ranging from $T_{min}=1.53$ to $T_{max}=3.28$. This naturally translates to $92\times 92=8464$ samples for the stacked Ising model at each temperature pair $(T,\tilde{T})$.

For the Rydberg system, we make use of a recent Quantum Monte Carlo method [76] to generate occupation basis samples of the Rydberg Hamiltonian. The QMC simulation is based on a power-iteration scheme which projects out the ground-state of the Hamiltonian:

$$\ket{E_{0}}\approx(C-H)^{M}\ket{+}^{\otimes N},$$

(4)

where $\ket{+}$ is the positive eigenstate of $\sigma_{x}$, $N$ is the number of lattice sites, $C$ is a constant energy shift used to cure the sign-problem emerging from the diagonal part of the Hamiltonian, and $M$ is called the projection length. We perform our simulations on a $16\times 16$ square lattice with open boundaries at various parameter values. Unlike previous DMRG-based studies[75, 69] we do not impose a truncation on the vdW interaction. We take our projection length $M$ to be 100,000 which we found was more than enough to accurately converge to the ground-state over the parameter sets which were simulated. For our simulations we fixed $\Omega=1$ and performed scans over $R_{b}\in\{1.0,1.1,\ldots,1.8\}$ and $\delta\in\{0.5,0.6,\ldots,2.9\}$.

To generate the occupation basis data for the SNNs, we first perform 100,000 Monte Carlo update steps to allow the chain to reach equilibrium. We then record one sample every 10,000 steps in order to eliminate any possible autocorrelation between successive samples. Each Monte Carlo step consists of a diagonal update step, followed by a cluster update step in which all possible line-clusters are built deterministically and flipped independently according to a Metropolis condition; see [76] for further details. Additionally, at each point in parameter space $(R_{b},\delta)$ we run 3 independent Markov chains. The chains are allowed to evolve until each has generated 400 samples, giving a total of 1200 independent samples for each parameter pair.

Artificial neural networks are directed graphs that have the ability to learn an approximation to any smooth function $f(x)=y$ given sufficiently many parameters. A neural network is built by successively applying matrix multiplications characterized by weights $w_{ij}^{L}$ that are offset by biases $b_{i}^{L}$. ($i,j$ are neuron indices in different layers $L$). Between subsequent matrix multiplications there is a non-linear activation function, common choices of which are sigmoid or rectified linear units. A neural network is trained by applying it to a data set and optimizing the network parameters to minimize a certain objective function using gradient descent.

Siamese neural networks (SNN) were introduced to solve an infinite class classification problem as it occurs in finger print recognition or signature verification [64, 65]. Instead of assigning a class label to a data instance, the SNN compares two data points and determines their similarity. A solution to the infinite class classification problem is obtained by calculating the predicted similarity between labelled anchor data points and a new unlabelled data instance.

A SNN $F(S_{i},S_{j})$ (Fig. 3) consists of two identical sub-networks $f$ (Fig. 2) which project an input pair into a latent embedding space. The similarity of two inputs is determined based on the distance in embedding space.

$\displaystyle F(S_{i},S_{j})=d\left(f(S_{i}),f(S_{j})\right)$

(5)

Possible distance metrics $d$ must be chosen according to the problem at hand, in our case we chose the average squared Euclidean distance on the unit sphere which is equivalent to the cosine distance.

$$d(f_{1},f_{2})=\frac{\sum_{i=0}^{N}(f_{1i}-f_{2i})^{2}}{N},$$

(6)

Here $f_{1i},f_{2i}$ denotes different components in an $N$ dimensional embedding space. Instead of training the SNN on paired training data, an effective way to train Siamese neural networks is through minimizing contrastive loss functions involving triplets $(S_{a},S_{+},S_{-})$ of data points

$$\mathcal{L}=\max(d(f(S_{a}),f(S_{+}))-d(f(S_{a}),f(S_{-}))+\alpha,0),$$

(7)

where the hyperparameter $\alpha$ is chosen such that the neural network is prevented from learning trivial embeddings. The intuition behind $\alpha$ is its interpretation as a margin to encourage separation between the anchor and the negative embedding in the embedding space [78]. Minimizing $\mathcal{L}$ requires minimizing the distance between the anchor and positive sample $d(f(S_{a}),f(S_{+}))$, while maximizing the distance between the anchor and negative sample $d(f(S_{a}),f(S_{-})$. $\max(\cdot,0)$ prevents the neural network from learning runaway embeddings $f(S_{i})_{i}\rightarrow\infty$.

The explicit model architecture for $f$ Fig. 2 depends on the underlying data set. In the case of the Ising model, $f$ consists of two 2-D convolutional layers, both with stride $(1,1)$, a kernel size of $(3,3)$, and with 6 and 10 filters, respectively. Each layer is fed into a ReLU activation function. The resulting image dimensions are $(16,16)$. This is followed by a $(16,16)$ average pooling layer. The output is flattened and fed to a dense layer with 10 neurons. Subsequently, we feed this output to the embedding layer, which also contains 10 neurons and a sigmoid activation function. The embedding is normalized to unit length under euclidean norm.

The model architecture for examining the Rydberg system is similar to the Ising model. In this case, we begin with two 2-D convolutional layers, both with stride (1, 1), a kernel size of (3, 3), and with 6 and 4 filters, respectively. The resulting image dimensions are (12, 12). As before, we subsequently apply an average pooling layer, but this time of dimension (12, 12). The remainder of the architecture is identical to before.

We use the Adam optimizer to train our neural network. Furthermore, our training procedure involves the early stopping callback. This technique involves training as long the loss is decreasing. If the loss is not decreasing for $m$ epochs, training is stopped. The value of $m$ is known as the patience. We set the maximum number of epochs to 150, which is enough to allow the callback to decide when to terminate the training process. We observe that small values of patience, around $m=1..3$ and $m=6$ are sufficient for training on the Ising model and the Rydberg system, respectively. Additionally, we use a learning rate of $\alpha_{lr}=0.001$ for the Ising model, and $\alpha_{lr}=0.0005$ for the Rydberg system. Another hyperparameter is the margin for the contrastive loss, where we use $\alpha=0.4$. We employ the TensorFlow and Keras libraries to implement the network, training, and callbacks.

In this scheme, a configuration corresponding to a point in the phase diagram is compared to its neighbors within a certain distance. At the example of a single Ising model, this means a configuration from temperature $T_{n}$ is compared to a configuration at $T_{n+k}$, where $n$ is the temperature list index and $k$ is a constant shift. The value of $k$ can be treated as a hyperparameter. A comparison between configurations at $T_{n}$ and $T_{n+k}$ is assigned a temperature of $\frac{1}{2}(T_{n}+T_{n+k})$.

Our similarity measure of choice is the normalized cosine dissimilarity scaled to a range $[0,1]$, where 1 corresponds to the empirically found maximal dissimilarity and 0 to the empirically determined minimal similarity. The cosine dissimilarity is related to the cosine distance via $1-cos(\theta)$, where

$$cos(\theta)=\frac{f(S_{n})\cdot f(S_{n+k})}{\lVert f(S_{n})\rVert\lVert f(S_{n+k})\rVert}$$

(22)

The cosine distance is equivalent to the euclidean distance when acting on the SNN embedding space normalized to the unit sphere ($\lVert f(S_{i})\rVert^{2}=1$):

$\displaystyle\lVert f(S_{n})-f(S_{n+k})\rVert^{2}$

$\displaystyle=2-2\ cos(\theta)$

(23)

Since the neural network $f$ is designed to output positive values $f(S_{n})_{i}\in[0,1]$, the cosine dissimilarity takes on its minimum $1-cos(\theta)=0$ for similar embeddings and its maximum $1-cos(\theta)=1$ if the embeddings are dissimilar.

The result of scanning across a phase transition is depicted schematically in Fig. 4(c). The highest peak will indicate the SNN prediction of the phase transition.

We examine the results of applying SNN unsupervised phase recognition to the $20\times 20$ stacked Ising model. Analytically, the phase boundaries of the stacked Ising model in the thermodynamic limit are $(T,\tilde{T})=(2.269,\cdot)$ and $(T,\tilde{T})=(\cdot,2.269)$. The phase transition temperature is prone to finite size effects as the lattice becomes smaller. If the phase transition would be calculated using the magnetization, finite size effects would distort the phase boundaries to $(T,\tilde{T})\approx(2.36,\cdot)$ and $(T,\tilde{T})\approx(\cdot,2.36)$; however, it is important to note that the quantities a neural network learns might differ from the magnetization and thus experience different finite size scalings [55].

In order to calculate the stacked Ising model phase boundaries, we choose to perform vertical and horizontal scans across the phase diagram, where each scan is repeated five times and the standard deviation of this ensemble is displayed as uncertainty. Exemplarily, the results of two of these scans can be found in Fig. 5. The location of these scans is depicted in (c). (a) reveals the phase transition from (ferromagnetic, paramagnetic) to (paramagnetic, paramagnetic) (FP to PP), while (b) displays the phase transition from PF to PP. Collecting 21 vertical and 21 horizontal scans reveals the phase boundaries of the stacked Ising Model, as seen in Fig. 6. By comparing the horizontal (green dotted line) and vertical scans (red dotted line) with the Ising model phase boundaries in the thermodynamic limit (black dashed lines) one can observe a clear difference. The SNN predicts critical temperatures of $T_{c}\approx 2.352$, consistent with the finite size correction of the magnetization which indicates a phase transition at $T_{c}\approx 2.36$.

In order to reveal the power of our SNN based phase recognition scheme, we have to scan across diagonal slices within the phase diagram that contain more than one phase transition. For this purpose we scan diagonally across the phase diagram as depicted in Fig. 7. The two diagonal scans across the lattice reveal that this technique can identify more than one phase transition. The first diagonal scan is performed along a line through the center of the lattice. In this case, we see a single phase transition, as the network scans directly through the intersection of the vertical and horizontal phase lines. A closer examination of the effect of changing $k$ can be found in Fig. 20. The prediction yields a phase boundary at $(T,\tilde{T})_{c}=(2.28,2.28)$.

The second diagonal scan we perform is shifted, such that it crosses both the vertical and horizontal phase transition. In this case, the true phase boundaries are cut by the shifted diagonal slice at $(T,\tilde{T})_{c}=(1.827,2.269)$ and $(T,\tilde{T})_{c}=(2.269,2.711)$. The network is able to capture both transitions at $(T,\tilde{T})_{c}\approx(1.857,2.299)$ and $(T,\tilde{T})_{c}\approx(2.281,2.723)$. A closer examination is found in Fig. 21.

Furthermore, the embeddings shed some light on how the neural network encodes the aforesaid phase transitions. In Fig. 8 we see embeddings which encode phase information for the diagonal slice. These embeddings clearly separate between the two underlying phases. Embedding 1 only spikes in the PP phase, while embedding 10 activates at FF regimes. A switch between both embeddings occurs at the phase transition. Other embeddings can be found in Fig. 16. In the case of a single embedding, we may have noise contained in the same embedding neuron where the phase behaviour is captured. With additional embeddings, the noise may be relegated to other embedding neurons, isolating the phase transition information.

The embedding space of the shifted diagonal scan is depicted in Fig. 9. We observe three crucial behaviours in the embeddings. Each of the three histograms show maximum activity in three different regions. Embedding 2 encodes the FF phase, embedding 4 the PP phase, and embedding 7 encodes configurations from the FP phase.

The Rydberg atom array phase diagram is not as well studied as the Ising model, so we need to compare our SNN phase boundaries to recent papers [69, 75] and evaluate the order parameters on the QMC data ourselves.

In order to identify the approximate phase boundaries we compute predictors for the phases of interest. Each phase corresponds to various peaks in the absolute value of the Fourier transform (FT) of the one-point function:

$$n(\mathbf{k})=n(k_{x},k_{y})=\absolutevalue{\frac{1}{\sqrt{N}}\sum_{j}n_{j}\exp(i\mathbf{k}\cdot\mathbf{r}_{j})}$$

(24)

where $n_{j}$ is the Rydberg occupation of the $j$th site. Furthermore, we symmetrize the FT by averaging over permutations of the momentum axes:

$$\mathcal{F}(k_{x},k_{y})=\frac{1}{2}\left[n(k_{x},k_{y})+n(k_{y},k_{x})\right]$$

(25)

The peaks occur for the checkerboard phase at $(\pi,\pi)$, for the striated phase at $(\pi,0)$ and $(\pi,\pi)$, and for the star phase at $(\pi,0)$, and $(\pi/2,\pi)$[75]. We compute, for each point in parameter space $(R_{b},\delta)$, the symmetrized FT at these specific momenta averaged over the full 1200 sample data set given to the SNNs.