arXiv CSopen access

TetrisCNN for interpretable detection of phases of matter from experimental quantum simulator data

· arxiv_cs

arXiv CS · Papers · License: Open Access

machine learning, deep learning, neural networks

TetrisCNN for interpretable detection of phases of matter from experimental quantum simulator data Kacper Cybiński,1, 2, ∗ Björn van Zwol,3, ∗ James Enouen,4 Guillaume Bornet,5 Thierry Lahaye,6 Antoine Browaeys,6 Antoine Georges,7, 8, 9, 10 and Anna Dawid3, †

arXiv:2609.20693v1 [quant-ph] 17 Sep 2026

Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland 2 IDEAS Research Institute, Królewska 27, 00-060 Warsaw, Poland 3 Applied Quantum Algorithms ⟨aQaL ⟩, LIACS & LION, Leiden University, The Netherlands 4 Department of Computer Science, University of Southern California, Los Angeles, CA 90089, USA 5 Princeton University, Department of Electrical and Computer Engineering, Princeton, New Jersey 08544, USA 6 Université Paris-Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire Charles Fabry, 91127 Palaiseau Cedex, France 7 Collège de France, PSL University, 11 place Marcelin Berthelot, 75005 Paris, France 8 Center for Computational Quantum Physics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA 9 CPHT, CNRS, École Polytechnique, IP Paris, F-91128 Palaiseau, France 10 DQMP, Université de Genève, 24 quai Ernest Ansermet, CH-1211 Genève, Switzerland (Dated: September 18, 2026) Detecting phases of matter in general relies on identifying the correct order parameter - a task that remains notoriously difficult for unknown transitions and traditionally is guided by physical intuition and educated guess. Neural networks have recently offered an alternative route by locating phase transitions in known models without any a priori physical knowledge. Yet these approaches remain black boxes and only identify phases without elucidating their properties. Moreover, they often struggle when confronted with realistic, noisy experimental data, which constitute the ultimate testbed for automated methods in physics. Here, we bridge these perspectives by introducing TetrisCNN, a convolutional architecture with parallel branches of differently shaped filters, reminiscent of Tetris blocks, that learns sparse, interpretable latent representations directly in terms of spin correlators. Applied to experimental snapshots of two-dimensional Ising and XY quantum simulators measured in multiple bases, the network not only detects phase transitions and crossovers but also expresses its latent representation and decision boundaries as symbolic formulas built from experimentally measurable spin correlators. This framework opens the way to integrating interpretable neural networks with quantum simulators to uncover and understand new phases of matter.

INTRODUCTION

Understanding complex quantum systems with longrange interactions and strong entanglement remains a central challenge of modern quantum physics, particularly as classical simulation becomes computationally prohibitive. Programmable quantum simulators [1–4] provide a remarkable experimental alternative, granting access to regimes beyond the reach of conventional numerical methods and potentially uncovering new phases of matter thanks to their tunability [5]. However, identifying suitable order parameters for previously unknown phase transitions remains nontrivial and often depends on expert intuition. The task requires navigating an exponentially large Hilbert space and analyzing the symmetries of the system, guided by physical insight and educated guesses. In parallel, machine learning has an increasing impact on quantum sciences [6–10], mirroring its influence on industry [11]. In particular, neural networks have emerged as powerful tools for detecting phases of matter [12–20], and they are especially promising in unsupervised settings, where no prior knowledge of the system or its order

∗ These two authors contributed equally. † [email protected]

parameters is assumed. In principle, such techniques could enable the direct identification of novel and exotic phases from experimental data, including ones that may have escaped human intuition. In practice, however, machine learning has uncovered few unexpected phases, with notable examples using autoencoders [21] and support vector machines [22]. Such discoveries remain rare partly because machine learning is seldom applied directly to experimental quantum-simulator data [10]. One reason for the challenging nature of learning from experiments is that, in practice, experimental constraints often restrict the set of accessible observables [23], preventing access to full information about the quantum state. Moreover, real-world measurements introduce noise and finite-size effects that can obscure phase transitions, and put limits on system sizes [24]. These challenges complicate both training and the interpretation of results, and have made supervised learning the most common approach when analyzing experimental data [25–29]. However, supervised models rely on ground truth and struggle with out-of-distribution generalization, limiting their potential for scientific discovery. Additional concerns include the stability of learning models against noise in experimental settings that can act like adversarial perturbations [30] and therefore fatally derail network predictions [31]. Even when unsupervised detection schemes are success-

FIG. 1: TetrisCNN processes snapshots through parallel branches with differently shaped filters that map to spin correlators. Sparsity suppresses redundant or task-irrelevant branches, yielding a low-dimensional interpretable bottleneck from which the decision boundary can be expressed symbolically in terms of spin correlators.

fully applied directly to experimental data [24, 27, 29, 32– 35], the resulting machine learning models typically provide little physical insight, with only a few notable exceptions [29, 33–35]. Accurate prediction alone does not encourage neural networks to find interpretable descriptions: they lack an intrinsic drive to “compress knowledge” [36–39]. Extracting reliable physical insight is further complicated in experimental settings, where explainability methods can be even more sensitive to noise than the underlying predictions [40, 41], while the Rashomon effect (the coexistence of multiple, equally predictive yet qualitatively distinct models, whose number grows with noise) undermines the stability of inferred conclusions [42]. An ultimate goal of machine learning for phases of matter is therefore to make such automated approaches interpretable [41, 43–51], especially for experimental data, enabling not only the detection of phases of matter but also an understanding of their physical properties, in particular the underlying order parameters [29, 35, 52–63]. To address these challenges, we introduce TetrisCNN, an interpretable convolutional neural-network architecture for identifying the spin correlators underlying phase transitions and crossovers directly from experimental snapshots. The key idea is to provide the network with an explicit drive to compress knowledge. TetrisCNN uses parallel convolutional branches that represent candidate local correlators at different spatial scales and geometries, with the resulting collection of differently shaped filters motivating the name “TetrisCNN”. Sparse regularization then automatically selects a small subset of correlators sufficient for a given learning task. This design enables the latent representation and the resulting decision boundaries to be expressed analytically in terms of experimentally measurable spin correlators, providing a missing ingredient of neural network-based approaches: interpretability. We apply TetrisCNN to experimental two-dimensional (2D) Rydberg-array data that realize the Ising [64] and XY [65] models, where it recovers sparse, physically mean-

ingful representations of the observed transitions and reveals the correlators underlying the network predictions. Throughout this work, we use transition broadly for pronounced changes in the system’s behavior, encompassing both phase transitions and crossovers. Beyond accurately detecting transitions, TetrisCNN reveals why the predictions are made, allowing machine-learning conclusions to be analyzed and challenged using theoretical arguments. The manuscript is organized as follows. We first introduce the experimental Ising and XY datasets, the learning task, and the TetrisCNN architecture, including the sparse bottleneck and Boolean-Fourier mapping to spin correlators (Sec. II A–II C). We then present the results in Sec. III, discuss their physical interpretation in Sec. IV, and compare TetrisCNN with related interpretable approaches in Sec. V. We conclude in Sec. VI. Implementation details and all architectural and numerical settings are provided in Apps. A–I, together with the open-source code and experimental datasets at GitHub [66].

II. A.

METHODS Quantum data

We showcase the usefulness of TetrisCNN on the example of interpretable detection of phases of matter from raw experimental projective measurements (snapshots) taken on the Rydberg quantum simulator of two quantum spin Hamiltonians: 2D Ising and XY models. a. Rydberg quantum simulators Among the most promising quantum simulation platforms are Rydbergatom arrays [3, 67]. Atoms trapped in optical tweezers form highly controllable pseudospin- 12 qubits that can be individually positioned, addressed, and measured with single-site resolution. Quantum correlations and entanglement are generated via strong electric dipole–dipole interactions, giving rise to the Rydberg blockade mecha-

FIG. 2: Experimental quantum data. (a),(c) Experimental protocols for the Ising and XY datasets, including control parameters, measurement points, and representative snapshots from the 8 × 8 Ising [64] and 6 × 7 XY [65] Rydberg-atom quantum simulators. (b),(d) Selected snapshot-averaged spin correlators, cf. Eq. (3), with shaded regions indicating their standard deviation across the sweeps: mstag = CA [■] − CB [■], with A, B being the two X 2 sublattices, CNN = Crot [■■] = 12 (C[■■] + C[■ ■]), (m ) is in-plane ferromagnetic magnetization squared, cf. Eq. (15). nism that suppresses simultaneous excitation of nearby atoms. Crucially, the geometry, interaction range, and driving protocols are highly tunable and programmable, which (after accounting for experimental imperfections [68]) makes them exciting platforms for exploring emergent many-body phenomena and realizing and studying Hamiltonians that are inaccessible to classical numerical methods. b. 2D Ising Hamiltonian A 2D array of Rydberg atoms coupled via repulsive van der Waals interactions naturally realizes the Hamiltonian of the transverse-field Ising model (TFIM) when driven to Rydberg states [64]: ĤIsing =

X i<j

Uij ni nj +

X ℏΩ X X σi − ℏδ ni , 2 i i

(1)

where the Rydberg and ground states are the (pseudo)spin states |↑⟩ = |75S1/2 , mJ = 1/2⟩ and |↓⟩ = |5S1/2 , F = 2, mF = 2⟩, respectively, with F and mF denoting the hyperfine quantum numbers and mJ the magnetic quantum number of the fine structure. The interaction strength between atoms is given by the van der Waals potential 6 Uij = C6 /rij , where C6 is the van der Waals coefficient and rij is the distance between atoms i and j. Here ni = |↑⟩⟨↑|i = (1 + σiZ )/2, σi are the Pauli matrices, and ℏ is the reduced Planck constant. The two spin states are coherently coupled by a laser field with Rabi frequency Ω and detuning δ, which play the roles of the transverse

and longitudinal fields, respectively. For a square lattice with lattice spacing a = 10 µm and atoms excited to the n = 75 Rydberg state, the nearest-neighbor interaction strength is U/h ≈ 1.95 MHz.

For this Hamiltonian, snapshots are measured in the Z basis for t ∈ [0.6, 6] µs at intervals of 0.2 µs. During this protocol, δ and Ω are varied quasi-adiabatically, as shown in Fig. 2(a) (see Ref. [64] for details). The system is initially prepared in the Z-polarized paramagnetic ground state | ↓↓ . . . ↓⟩. As the detuning is increased while the transverse field remains finite, the system first undergoes a smooth crossover out of this nearly fully polarized state and enters a more strongly mixed paramagnetic regime. At later times, the transverse field is reduced and the protocol crosses the paramagnetic-to-antiferromagnetic quantum phase transition, leading to the growth of staggered magnetization.

c. XY Hamiltonian Another experimental setup [65], whose measurements we study, consists of a 2D rectangular lattice of 87 Rb atoms in an optical 6 × 7 tweezer array, where an effective spin-1/2 is encoded in the Rydberg states |↑⟩ = |60S1/2 ⟩ and |↓⟩ = |60P1/2 ⟩. Resonant dipole– dipole interactions allow the realization of a long-range dipolar XY Hamiltonian supplemented by a staggered

4 longitudinal field, leading to the total Hamiltonian ĤXY,δXY = −

J X a3 X X σi σj + σiY σjY 3 2 i<j rij

+ ℏδXY

X 1 − (−1)ix +iy i

(2)

ni .

Here, J/h = 0.77 MHz, the lattice spacing is a = 12.5 µm, and a magnetic field perpendicular to the lattice plane defines the quantization axis, resulting in isotropic dipolar interactions. The integers (ix , iy ) denote the coordinates of site i on the square lattice, such that (−1)ix +iy alternates between the two sublattices. With the convention (−1)ix +iy = +1 on sublattice A and −1 on sublattice B, the staggered-field term vanishes on A and equals ℏδXY ni on B, where ni = (1 + σiZ )/2. It therefore describes the staggered light shift used experimentally to prepare the initial Néel state. For a sufficiently large positive value of δXY , the initial Néel state approximates the ground state of ĤXY,δXY . Reducing the staggered field toward zero following an approximately adiabatic exponential ramp, as shown in Fig. 2(c), connects this state to a low-energy state of the XY Hamiltonian with ferromagnetic correlations in the XY plane. Conversely, for a sufficiently large negative staggered field, the same Néel configuration approximates the highest-energy state, and the corresponding ramp produces antiferromagnetic XY correlations [65]. In this work, we analyze the ferromagnetic protocol. Snapshots are measured in both the X and Z bases for t ∈ [0, 8] µs at irregular intervals. d. Working with experimental data vs simulated data Most deep learning studies in quantum physics rely on numerically simulated ground states. TetrisCNN was also initially developed as a proof of concept on numerically simulated datasets, including the one-dimensional (1D) TFIM and the 2D Ising lattice gauge theory [69]. We believe, however, that applying machine learning methods to experimental or noisy data is the ultimate test for their utility and robustness [41]. Therefore, here we analyze raw snapshots from quantum simulation experiments. These configurations are produced by finite-time approximately adiabatic sweeps and therefore represent non-equilibrium dynamical states rather than exact ground states. The data additionally contain realistic experimental imperfections, including state preparation, detection, and basisrotation errors, inhomogeneous fields, interaction disorder from atomic-position fluctuations, and decoherence during the ramps. While local errors lead to bounded errors in local observables, these imperfections reduce the fidelity of the quantum many-body state and pose increasing challenges for preparing and characterizing larger systems. Since the experiments are performed on finite lattices with open boundaries, edge and finite-size effects also influence the observed ordering patterns. Finally, comparisons with classical equilibrium ensembles [64] reveal systematic biases toward more ordered configurations, indicating that

the data reflect specific noisy quantum dynamics rather than thermal or idealized equilibrium samples. e. Datasets We now take a machine learning perspective. Projective measurement snapshots from the models above are taken as the network input x(n) = {Sib }D i=1 , with i indexing the spatial grid and b ∈ {X, Z} denoting the measurement basis. When both bases are available, we randomly combine snapshots acquired at the same sweep point into two channels as x(n) = {(SiX , SiZ )}D i=1 . Each snapshot has a corresponding label y (n) that depends on the learning task, which we specify further in Sec. II B. Snapshots and labels are collected in a dataset D ≡ {(x(n) , y (n) )}N n=1 . We focus on 2D square grids, so the data is represented as a tensor of shape [N, C, D1 , D2 ] with D = D1 × D2 the number of spatial sites, C = 1 by default, and C = 2 if we consider two bases simultaneously. D is divided into training and validation sets in a 7:3 ratio. We measure evaluation metrics on the validation set for early stopping, and to compare the fits of different models. We do not use a test set because we are not interested in benchmarking; rather, we use machine learning as a supervised subroutine in a broader unsupervised discovery process (i.e., evaluation metrics on the validation set should not be interpreted as expected model performance, as is common in machine learning). More information on dataset preparation is in App. A. f. Order parameters In general, quantum systems exhibit a wide variety of types of order. The Landau paradigm defines order by symmetry breaking [70, 71] (e.g. in the TFI and XY model above); BKT order relates to binding/unbinding of topological defects [72, 73] (e.g. the AFM variant of the XY model); topological order is characterized by long-range entanglement and topologydependent ground state degeneracy [74, 75] (e.g. the toric code [76]). Order can be even more unconventional [77, 78] and novel types of order may yet be discovered. We desire a machine learning approach that provides insight regardless of the system’s order – here, we take a step towards this goal. g. Spin correlators For a machine learning algorithm to provide insight, it needs to be expressed in a language we can recognize. A general way to characterize order in many-body systems is through spin correlators or npoint functions. For example, ferromagnetic order can be characterized by the one-point function, or magnetization, Si = m, where the overline bar denotes a spatial average over indices i – or using long-range behavior of the two-point function lim|i−j|→∞ Si Sj . Staggered magnetization can be written as SA − SB where A and B are two sublattices. In general, correlators are characterized by the number of spins they involve. More complicated phases, however, involve correlators in specific patterns (e.g. staggered, rhombic or striated patterns [79, 80]). We therefore represent a spin correlator by a corresponding geometrical pattern P that specifies the relative position of spins whose product is taken. For a given snapshot, let TP denote all translations of P that fit

5 within the physical system. We define the corresponding single-snapshot spin correlator as C[P ](S) =

1 X Y Si |TP | T

(3)

n,i

i∈T (P )

that is, the product of the spins within a pattern P translated by T . Throughout, sums over T are always defined over all possible translations of the pattern within the snapshot, T ∈ TP . From here onwards, we shall omit the dependence on spins (S) to avoid clutter. For example, P = ■■ = {0, 1} defines the horizontal nearest-neighbor spin pair, and therefore P 1 S S C[■■] = D−1 = S S i i+1 , for a snapshot of a i i i+1 one-dimensional system of size D. Similarly, the singlesite pattern P = ■ gives the snapshot magnetization, C[■] = m. When no lattice orientation is physically distinguished, it is useful to consider correlators that are invariant under lattice rotations. We define the rotationally invariant correlator associated with P by averaging over its distinct rotations on a square lattice, 1 X Crot [P ] = C[R], |RP |

(4)

where R ∈ RP are all distinct 90◦ rotations of P , corresponding to the C4 symmetry group. For example, Crot [■■] = 21 C[■■] + C ■ ■ . Throughout this work, C[P ] denotes a quantity evaluated on a single snapshot. We denote its expectation value over experimental snapshots by ⟨C[P ]⟩ (and analogously m and ⟨m⟩ for the magnetization). Several such correlators are plotted in Fig. 2(b) and (d) for the TFIM and XY datasets. B.

When using TetrisCNN for supervised learning, we employ the cross-entropy loss: X (n) (n) LCEL = − yi log ŷi ,

Task

Several unsupervised methods for finding phase transitions have been proposed in the literature. These include clustering techniques applied within a low-dimensional space, which can be obtained by dimensionality-reduction techniques [81], diffusion maps [82], or autoencoders [14, 21, 33, 34, 83, 84]. Other approaches include learning by confusion [13], prediction-divergence method [15, 46, 85], generative modeling [18], and distance learning [20, 24]. TetrisCNN is an architecture that can serve as a drop-in replacement in a large number of methods: autoencoders, prediction divergence (App. F 6), learning by confusion, distance learning, and supervised classification. For clarity and simplicity, we present in this work how TetrisCNN can be used for supervised classification. The estimated phase boundary is determined separately in a prior stage: we use and compare several methods and adopt the transition point that is most frequently observed among these predictions.

with i ∈ {0, 1} the class label and ŷ a softmax of the final layer. Task performance is measured by accuracy, the fraction of snapshots correctly classified. Note that the network here makes predictions on single snapshots, which sometimes limits the maximum accuracy it can achieve. For example, if the same spin configuration appears in different phases, the network cannot make correct predictions for all of them. Moreover, networks trained on single snapshots cannot compute directly so-called ensemble connected correlators [86] such as ⟨Si Si+1 ⟩c = ⟨Si Si+1 ⟩ − ⟨Si ⟩⟨Si+1 ⟩, as in the in-plane magnetization squared (mX )2c in Fig. 2(d), cf. Eq. (15). C.

Interpretable architecture of TetrisCNN

a. Interpretability and expressivity Neural networks (NNs) stand out among machine learning methods by being both trainable and highly expressive. This expressivity is associated with network depth, width, and total parameter count [87, 88]. Neural networks also excel at representation, or feature learning [89, 90]. Extracting or interpreting these features from the network, however, is typically extremely difficult. Hence, it is commonly assumed that expressivity and interpretability are inherently opposed [91–93]. Our network design, TetrisCNN, showcases how this is not necessarily true for Boolean-valued data (e.g., spins). This is achieved using a combination of elements that conspire to enable interpretation – here defined in a broad sense. We briefly explain each design element below, including its contribution to interpretability. We emphasize the generality of our approach; TetrisCNN provides a drop-in replacement for any neural network that was heretofore a black box. b. Locality and spatial symmetries Our data is defined on a 2D square lattice, implying an approximate translational symmetry (an infinite lattice would be fully symmetric). The physics defined by Eqs. (1) and (2) is 6 3 also mostly local, interactions decaying as 1/rij and 1/rij , respectively. This motivates the use of convolutional layers, which enable highly efficient learning on data with these properties [94]. Reusing the notation in Eq. (3), we express a convolution in a non-standard form using a pattern P : X conv[P ]T = Wj ST (j) . (5) j∈P

Here, S is the full input spin configuration (snapshot), while pattern P specifies the support of the convolution. We again omit the dependence on S on the left-hand side to avoid clutter. Wj is its associated trainable weight,

1x1 1x2 0

0x0 1x1 1

1x1 0x2 0

1x0 1x1 1

2.67

Global Average Pooling (GAP)

...

2D input data

Convolution

Feature map

Filter 2x2

FIG. 3: Example of convolution with a 2 × 2 filter, followed by global average pooling (GAP). The filter is applied across the input to produce a feature map, which GAP reduces to a scalar.

and T (j) is the corresponding translated lattice site. The pattern P , together with the weights {Wj }j∈P , defines a filter : a local function that is applied with the same weights at every allowed translation T , as illustrated in Fig. 3 for a 2 × 2 filter. Thus, Eq. (5) is the filter output at translation T and depends only on the spins, and {Si }i∈T (P ) . Applying the same filter at every allowed translation of P within the snapshot produces a feature map. Convolutions are typically paired with pooling layers, which aggregate or downsample the resulting feature map. Of particular interest here is global average pooling (GAP), which averages the feature map over all spatial positions (see the final arrow in Fig. 3). The convolution in Eq. (5) is equivariant: translating the input translates the feature map. GAP then removes the spatial index, making the pooled output translation invariant (up to finite-size boundary effects). While standard convolutions (Eq. (5)) capture translational symmetry, the C4 rotational symmetry of the square lattice can be accounted for with group-convolutions. Following [95], we implement this by applying each learned filter in all distinct rotations R ∈ RP , rotating its pattern and weights together, and averaging the resulting feature maps: convrot [P ]T ≡

1 X conv[R]T . |RP |

(6)

Here, R ∈ RP denotes all distinct 90◦ rotations of P , and conv[R] denotes the correspondingly rotated version of the same filter, including its weights. This makes the convolutional feature map rotationally equivariant, and its spatially pooled output rotationally invariant. c. Interpreting convolutions with correlators Previous work [52] showed that when any function of a convolution (e.g. a convolution followed by a nonlinearity) is combined with global average pooling (GAP), a remarkable interpretation can be made for Boolean-valued spin

configurations S ∈ {±1}D Let f [P ]T denote the output at the translated pattern T of such a local function with receptive field P . Using a Boolean Fourier expansion [96–99] (cf. App. B 1), one can show that X 1 X cP ′ (f ) C[P ′ ] , f [P ]T = (7) |T | P

P ′ ⊆P

where cP ′ (f ) ∈ R are coefficients determined by the function f [P ], and P ′ ⊆ P runs over all subpatterns of P . Thus, any local function of Boolean-valued spins, when combined with GAP, produces a linear combination of the spin correlators C[P ′ ] defined in Eq. (3). This was previously shown in Ref. [52] for small patterns (■ and ■■), e.g., for P = ■■, the GAP output takes the form c∅ + c1 C[■] + c2 C[■■] (up to open boundary effects). Eq. (7) generalizes this result to arbitrary patterns (correcting a minor theoretical error in their derivation in the process, replacing it by a principled foundation using Boolean Fourier analysis). It should be noted that {cP ′ } are not directly available, but can be obtained post hoc by simply taking a linear regression fit on the output using {C[P ′ ]}. This decomposes the learned function into physically meaningful spin correlators. We also note that this procedure is efficient up to a 3 × 3 filter but becomes expensive for larger filters, due to exponential scaling of the number of subpatterns P ′ ⊆ P as |P | grows (400 and 57,856 translationally symmetric subpatterns for 3 × 3 and 4 × 4 filters, respectively; see App. B 2). d. Interpretable branch definition With this result, we now define a TetrisCNN branch as i 1 Xh z[P ] ≡ conv[■] ◦ ϕ ◦ conv[P ] . (8) |TP | T T

This sequence of operations is designed to maximize expressivity through layers and channels, while preserving the expansion in Eq. (7). Specifically, the first convolution conv[P ] sets the branch’s receptive field to P , and produces 32 output channels. The subsequent 1 × 1 convolution, conv[■], combines these channels locally without enlarging the receptive field. The function ϕ is an elementwise nonlinearity, chosen here to be ReLU. Finally, GAP averages the resulting feature map over all translations T , producing the scalar branch output z[P ]. As examples, we show several z[P ] in Tab. I for patterns of increasing size, along with the correlators that each z[P ] can compute. e. Low-dimensional latent space A common approach for obtaining interpretable features is to project data onto a latent space before making a prediction. By further enforcing it to have specific properties (e.g., a given dimensionality or sparsity), one provides an explicit drive to compress knowledge, and allows one to ‘disentangle’ the neural network’s representation [100–102]. We employ this strategy in TetrisCNN, although in a non-standard way. That is, we use a composition of mappings: TetrisCNN

Task NN

x 7−−−−−−−−→ {zk }k=1 7−−−−−−→ ŷ, | {z } Bottleneck

(9)

7 Correlators C[P ′ ], P ′ ⊆ P C[∅] = const. z1 conv[■]i = W0 Si C[■] = Si C[∅] = const. C[■] = Si z2 conv[■■]i = W0 Si + W1 Si+1 C[■■] = Si Si+1 C[∅] = const. C[■] = Si z3 conv[■■■]i = W0 Si + W1 Si+1 + W2 Si+2 C[■■] = Si Si+1 C[■□■] = Si Si+2 C[■■■] = Si Si+1 Si+2 Convolution filter conv[P ]

TABLE I: Mapping filters to spin correlators. A TetrisCNN branch zk computes a linear combination of all correlators C[P ′ ] for subpatterns P ′ ⊆ P , defined from a convolution conv[P ] (which follows from a Boolean Fourier expansion and GAP, see Sec. II C and App. B 1). By additionally promoting Psparsity on branch outputs using L1-regularization k λk |zk |, where λ1 < λ2 < ... < λK , TetrisCNN selects minimally sized task-relevant correlators. E.g. if trained with z1,2,3 above and z3 is the only non-vanishing branch output, {∅, ■, ■■} are task-irrelevant, since they are included also in z2 which vanished during training despite weaker penalization, λ2 < λ3 . Thus one is left with ■□■ and/or ■■■. Note that an unfilled square denotes a masked-out site. For rotation-averaged convolutions convrot , the corresponding spin correlators Crot are likewise rotationally invariant. K

where {zk }k=1 is the latent space, also called bottleneck. In typical autoencoder architectures, x → 7− z is a fully connected neural network. In TetrisCNN, it instead consists of K parallel branches x → 7− zk , with each branch zk = z[Pk ] defined by Eq. (8). This design achieves a second level of interpretability compared to standard autoencoders. Namely, such bottlenecks are only ‘interpretable’ in the sense that they have low dimensionality – each individual dimension is still a complicated (non-interpretable) function of the data. By contrast, TetrisCNN additionally achieves interpretability of individual bottleneck dimensions by virtue of Eq. (7). f. Sparsity regularization The chosen patterns control the correlations that TetrisCNN can learn. But, as visualized in Tab. I, including branches with larger filters creates redundancy in the correlations picked up by different branches. We break this redundancy by additionally making the bottleneck sparse through L1-regularization: LL1 =

λk |zk |,

(10)

where λk are real-valued hyperparameters. We also call this the branch penalty. Adding this term to the loss means the network suppresses branches that are not relevant for the task during training. This achieves automatic feature selection: branches whose activations remain above the optimization-induced floor (see App. D) are retained as task-relevant, while unnecessary branches are suppressed.

g. Penalizing large patterns The Rashomon effect predicts that multiple models with equally good performance exist for a given dataset [42], which we also observe in practice (App. F 1). To break this redundancy, we introduce on top of Eq. (10) a penalty for ‘complexity’, here defined simply as the cardinality of the pattern |P | (which translates to the degree of the corresponding largest correlator, a canonical measure of complexity [97]). We interpolate λk linearly in log space between λmin and λmax (in the main text equal to 10−3 and 103 , respectively) according to |Pk |: α λk = λ1−α min λmax

(11)

|Pk |−1 where α ≡ max |Pk |−1 ∈ [0, 1]. This means that branches with larger patterns incur a higher cost and are thus suppressed unless they provide additional predictive power compared to branches with smaller patterns. In this way, since Eq. (10) makes the bottleneck sparse, TetrisCNN is biased toward the smallest pattern that provides predictive information sufficient to solve a task. We will refer to suppressed branches as ‘inactive’ or ‘deactivated’, and the remaining branches as ‘active’. In practice, we find that learning rate scheduling (App. D 3) helps accentuate the difference between active and inactive branches. h. Task network The innovative part of TetrisCNN relies on the interpretability of its bottleneck. These bottleneck activations are in the end input to a standard black-box fully connected network ({zk } 7− → ŷ), which we refer to as a task network (grey layer in Fig. 1). The degree to which the task network can be interpreted depends on the task. In the classification task, we treat the task network as a classifier operating in the bottleneck space, thereby yielding interpretable decision boundaries, as presented in the text below. In more complex settings, one can resort to symbolic regression (see App. I, also for criticisms of the method), which becomes cheap due to the bottleneck’s interpretability and sparsity. i. TetrisCNN hyperparameters TetrisCNN has two main hyperparameters. One is the maximum branch penalty λmax . Increasing λmax penalizes larger-filter branches more strongly and therefore favors solutions based on fewer and lower-order correlators. In App. F 1, we show that our results are robust to the specific choice of λmax , once a minimum threshold is reached. The second architectural hyperparameter is the choice of the set of patterns {Pk }. Here we take all subpatterns within a given 2 × 2 receptive field, Pk ⊆ ■■ ■■. We justify this in App. F 2, where another TetrisCNN run with larger filter sizes (1 × 1, 2 × 2, 4 × 4, 8 × 8) finds that scales beyond 2 × 2 are irrelevant for the tasks at hand. This procedure can also be viewed as a data-driven search for the spatial scale of the correlators relevant to the task. Moreover, we impose the C4 rotational symmetry of the square lattice, so each filter is applied in all distinct 90◦ rotations of its pattern, and the resulting feature maps are averaged. Consequently, each branch is invariant under rotations of its input pattern. In App. F 3 we show that rotationally invariant and unconstrained variants of TetrisCNN

8 ■

■■

■□ □■

(a)

■■ ■■

(b)

(c)

100

1.4

10−3

1.0

z[P ]

|z[P ]|

1.2

t (µs)

■■ ■□

10−6

−2

0.8 10−9 0

D L

P M A

A PC

PD M

0.6

1.0

10−1

0.8

10−3

(f)

D L

0.5

z[P ]

0.0 −0.5

A PC

M PD

(e)

10−9

0.2

t (µs)

10−7

0.4

150

10−5

0.6

100

Epoch

|z[P ]|

t (µs)

(d)

100

Epoch

200

−1.0

t (µs)

FIG. 4: Data-driven transition detection and interpretable classification from experimental snapshots for (a)–(c) Ising and (d)–(f) XY. (a,d) Transition locations estimated with established unsupervised methods (App. E). Red dashed lines mark the location used to define supervised labels. For Ising, the detected change is the polarization crossover rather than the later antiferromagnetic phase transition. (b,e) Evolution of rotationally invariant TetrisCNN branch activations during training, showing sparsity-induced suppression of task-irrelevant branches. (c,f) Activations of the surviving branches across the experimental sweeps. Unfilled squares denote masked-out sites.

give similar results. For further architectural details of TetrisCNN and information on the hyperparameters used, we refer the reader to App. C. The optimization details are in App. D.

III. A.

RESULTS

Sparse description of the identified transitions

a. Unsupervised detection of transitions across the experimental sweeps We first identify the locations of transitions in the two experimental datasets using several established unsupervised learning approaches, namely learning by confusion (LBC), the prediction-divergence method (PDM), principal component analysis (PCA), Uniform Manifold Approximation and Projection (UMAP), and distance learning (DL). Here and throughout, we use transition broadly to encompass both phase transitions and crossovers. The resulting location estimates are shown in Fig. 4(a,d) and discussed in detail in App. E. The transition location with the best agreement among

the methods is t = 1.1 µs for the Ising data, while the estimated transition location for the XY data is t = 0.625 µs. These locations are marked by red dashed lines in panels (a),(c)–(d), and (f) of Fig. 4. For all methods, the error bars span the acquisition times flanking the inferred transition. For LBC, they also reflect variation in the inferred transition location across runs. As discussed in Sec. IV, the change detected in the Ising data coincides with the loss of the initial Z-polarization rather than the later emergence of antiferromagnetic order. b. Supervised training of TetrisCNN Having identified the relevant transition in each dataset, we construct a supervised classification task by assigning different labels to snapshots from opposite sides of the identified location and train TetrisCNN on this task. The purpose of this step is not to improve the location estimate, but to reveal which local observables are most relevant for distinguishing the two classes. TetrisCNN reaches validation accuracies of 99.14% and 93.54% for the Ising and XY datasets, respectively. In both cases, most misclassifications occur near the identified transition (see Fig. 19(a),(b) in App. F 5), partly due to overlapping spin

Mapping branch activations to spin correlators

a. Boolean Fourier expansion A key advantage of TetrisCNN is that its interpretable architecture permits a direct interpretation of its latent representation in terms of physically meaningful spin correlators. To identify the observables encoded by the active branches, we express each branch activation as a multilinear expansion in the correlators associated with the corresponding filter pattern using the Boolean Fourier expansion, as described in Sec. II C; see Eq. (7). b. Magnetization in the Ising dataset For the Ising dataset, as shown in Fig. 4(b), only the z[■] branch remains active after training. From Tab. I, we see that the only correlators that a branch with a filter ■ can compute are C[■], so in this case Z-basis magnetization, C Z [■] = SiZ . A linear regression fit gives z[■] = −4.626 C Z [■] + 0.954 .

(12)

Thus, despite having access to a large collection of candidate correlators, TetrisCNN automatically identifies a 1D description of the phase classification problem. c. Three correlators are needed for XY data In the XY case, the network uses a richer set of correlators. As shown in Fig. 4(e), not only the z[■] branch but also the z[■■] and z[■□ □■] branches remain active, indicating that the task cannot be solved with high accuracy with a single local observable. In Fig. 5, the multilinear regression reveals that the z[■] branch is dominated by the X-basis magnetization C X [■], while the z[■■] and z[■□ □■] branches compute rotationally invariant Z-spin correlators, i.e., Z Z ■□ Crot [■■] and Crot [□■], respectively. The linear fits of each activation to its dominant spin correlator, given in full in Fig. 5(a)-(c), provide an almost exact description, with R2 = 1.00, 0.97, and 0.98, respectively.

z [] ≈ −1.56 C

(a) 0.5

X: Z:

+ 0.02 | R2 = 1.00

0.0 −0.5 −1.0 −1.5

X: Z:

(b)

X: Z:

z [] ≈ −0.21 Crot

0.1

X: Z:

∅

− 0.03 | R2 = 0.97

0.0 −0.1 −0.2 X: Z:

X: Z:

∅

≈ 0.16 Crot

X: Z:

(c)

X: Z:

other+ other−

− 0.01 | R2 = 0.98

0.2

0.1

configurations on opposite sides of the transition. c. Evolution of TetrisCNN branch activations Figures 4(b,e) show the evolution of branch activations during training. Due to the sparsity regularization, branches that are not required for the classification task gradually ‘deactivate’. For the Ising dataset, a single branch z[■] remains active after training, while all other branches are suppressed by several orders of magnitude. In contrast, multiple branches remain active for the XY dataset, indicating that a richer set of descriptors is required to distinguish the corresponding classes. The activations of the branches across the experimental sweep are shown in Fig. 4(c,f). In both datasets, the selected descriptors change sharply around the identified transition. In App. F 1, we show that the identified branches are robust to the choice of sparsity strength. Having established that only a small subset of branches is required for classification, we next identify the physical observables encoded by the active branches.

0.0

X: Z:

other+ other−

FIG. 5: Mapping from branch activations to correlators via multilinear regression for three active branches in TetrisCNN trained on the XY dataset as in Fig. 4(e)-(f), enabled by the Boolean Fourier expansion, Eq. (7). We need only to consider correlators C[P ′ ], P ′ ⊆ P defined by patterns P ′ that are contained in the filter pattern P (see Tab. I). Note that with all C[P ′ ], P ′ ⊆ P , the fit is perfect (R2 = 1). Small positive (negative) contributions are summed into “other+(–)”. The complete decomposition is shown in Fig. 21.

d. Variance of the z[■] branch Interestingly, the activation of the z[■] branch in Fig. 4(f) has much larger variance than other branches, particularly after the phase transition. We will resolve this riddle in the next section once we identify the function of the spin correlator used by the classifier through z[■]. Taken together, the Ising and XY results demonstrate that TetrisCNN not only identifies a sparse set of descrip-

−3

Decision boundaries in interpretable latent space

0.05z[■]2 − 0.03z[■] − 0.9z[■■] − z[■□ □■] + 0.03 = 0 . (13) Excitingly, we can use the mapping from activations to spin correlators in Fig. 5, and express this decision boundary in terms of spin correlators: 0.68C X [■]2 + 0.21C X [■] Z Z ■□ + Crot [■■] − 0.81Crot [□■] + 0.35 = 0 .

Phase 1 Phase 2

−2

−1

z []

(b)

0.2 0.1 0.0

−0.1

0.2

−2

0.0 −1

z [ ]

−0.1

]

0.1

a. Bottleneck space is interpretable So far, we have identified which spin correlators the network computes with its active branches to make the predictions. Now we can take a step further and understand how the network will classify a new snapshot by identifying the symbolic formula for its decision boundary. We achieve this by analyzing the network predictions directly in the latent space spanned by the branch activations. Crucially, this low-dimensional latent space is no longer an abstract embedding, as in, for example, autoencoders [14]: each axis has a known physical meaning in terms of spin correlators. b. The Ising bottleneck is 1D Disregarding the deactivated branches, we plot the activations of the only remaining active branch z[■] for all snapshots from the validation set in Fig. 6(a) and color code them by their predicted label. This reveals the full decision rule: the decision boundary is a single point at z[■] = −0.837, which, as we know from Eq. 12, corresponds to a Z-magnetization threshold of C Z [■] = SiZ = 0.387. c. The XY bottleneck is higher-dimensional Again, we plot all snapshots in the space of active branches, which for XY is three-dimensional (3D) and spanned by z[■], z[■■], and z[■□ □■]. Here, the decision boundary is a surface. We approximate it by fitting a polynomial surrogate to the TetrisCNN predictions, as detailed in App. H 1. We find that an almost perfect (99.84%) approximation is quadratic:

z [] = −0.837 C Z = 0.387

(a)

tors, but also reveals their physical interpretation in terms of experimentally measurable spin correlators.

FIG. 6: Decision boundaries in the interpretable latent space of TetrisCNN trained on the Ising (a) and XY (b) datasets. Each point represents one experimental snapshot, positioned according to its latent activations and colored by the phase label predicted by the classifier. (a) In the Ising case, only a single branch remains active, yielding a 1D latent space in which the classifier acts by thresholding the Z-basis magnetization. (b) In the XY case, the latent space is spanned by the three active branches identified in Fig. 4(e)-(f). The displayed surface is a quadratic approximation to the learned decision boundary, fitted on the training data, and reproduces the network’s predicted labels on the validation data with 99.84% accuracy. Because each latent coordinate corresponds to a known spin correlator (Fig. 5), both classification rules can be expressed directly in terms of physically meaningful observables, cf. Eq. (13) → Eq. (14).

(14)

which yields the final result: an explanation of the network’s classification rule in terms of physically relevant observables. The identified formula is robust across random initializations and random re-pairings of the X- and Z-basis snapshots, and we provide robustness tests and fitting details in App. F 4 and App. H 1, respectively. The identified quadratic contribution shows that identifying the relevant observables alone is not sufficient to explain the classifier. The task network learns a nontrivial function of these observables (such as z[■]2 ), which leads beyond a linear classification in the latent space. In this case, the quadratic contribution from C X [■] provides an improvement, especially at the phase transition where it reduces the disagreement with the TetrisCNN predictions from 3% to zero, as shown in App. H 2.

Supporting this interpretation, TetrisCNNs trained separately on the two measurement bases (App. F 5) recover a quadratic function of the X-magnetization from X-basis snapshots, while the Z-basis network recovers the same linear function of the two Z-spin correlators identified above. The latter provides substantially stronger singlesnapshot discrimination, reaching 91.19% validation accuracy compared with 70.28% for the X-only network. Thus, TetrisCNN combines the stronger overall singlesnapshot classification signal from the Z-basis correlators with the physically expected squared X magnetization, whose contribution is particularly important near the identified phase transition. d. Order parameter and the branch variance Remarkably, this quadratic term C X [■]2 closely relates to the

11 order parameter of the phase transition in the XY dataset identified by Ref. [65]; the in-plane (ferromagnetic) magnetization squared: (mX )2c =

1 X eX C N 2 i,j i,j

(15)

e X = σ X σ X − σ X σ X (not to be confused where C i,j i j i j with C[P ] from Eq. (3)). Rewriting, we obtain: (mX )2c = C X [■]2 − C X [■]

(16)

where TetrisCNN bases its decision boundary on the single-snapshot quantity C X [■]2 , whose average over snapshots gives the first term of the connected order parameter, C X [■]2 . Finally, we can also explain the variance of the z[■] branch in the XY dataset. The shaded region in Fig. 4(f) is: rD E p 2 (z − ⟨z⟩) = ⟨z 2 ⟩ − ⟨z⟩2 . Here, z = z[■]. Using the fit in Fig. 5(a), z[■] ≈ C X [■]. Thus, the squared width of the z[■] activation distribution is approximately the variance of C X [■], which is directly related to the in-plane ferromagnetic order parameter (mX )2 defined in Eq. (15). IV.

DISCUSSION

a. Polarization crossover in the Ising data All machine-learning approaches considered in this work, including learning by confusion, prediction divergence, principal component analysis, UMAP, and distance learning, consistently identify a transition around t ≈ 1.1 µs. The interpretability of TetrisCNN reveals that the corresponding classification signal is carried by the uniform Z-magnetization. This signal does not correspond to the later onset of antiferromagnetic order. Instead, it reflects the smooth crossover that occurs as the increasing detuning δ drives the system out of the nearly fully Z-polarized state and into a more strongly mixed paramagnetic regime while the transverse field Ω remains finite. b. Antiferromagnetic transition in the Ising data The later paramagnetic-to-antiferromagnetic quantum phase transition is signaled by the growth of the staggered magnetization in Fig. 2(b). On the finite experimental arrays, the transition is rounded, and the staggered magnetization increases continuously rather than displaying the sharp nonanalytic behavior expected in the thermodynamic limit. The system-size comparison reported in Ref. [64] shows that this increase becomes sharper for larger arrays, consistent with finite-size rounding, while the finite ramp duration and experimental imperfections introduce additional broadening. For reference, ground-state density-matrixrenormalization-group calculations in Ref. [64] locate

the inflection point of the staggered magnetization near t ≈ 3.2 µs for a 6 × 6 system and t ≈ 2.6 µs for a 10 × 10 system. Interestingly, none of the automated methods considered here detects this later transition, even though its signatures are visible in several correlators, Z including Crot [■■]. In contrast, a recent graph-theoretic data-driven approach reports a transition location consistent with the emergence of antiferromagnetic order [103]. Understanding why this broad class of methods preferentially detects the earlier polarization crossover and whether this preference is related to the different sharpness of the two changes remains an open question. c. The outcome of interpretation - or what we can learn More generally, the observables identified by TetrisCNN should be interpreted as those the network uses to solve a given learning task, rather than as straightforward order parameters or observables most fundamental to the physics of interest. The XY dataset provides an instructive example. Here, TetrisCNN identifies the quadratic contribution of the X-basis magnetization to the decision boundary, in close agreement with the established order parameter of Ref. [65]. However, this quantity alone does not provide sufficient information to reliably classify individual snapshots: a TetrisCNN trained only on X-basis snapshots reaches substantially lower accuracy than one trained on Z-basis snapshots (App. F 5). The network therefore supplements the physically expected squared X magnetization with two less conventional Z-basis correlation functions, C Z [■■] and C Z [■□ □■], which provide a much stronger single-snapshot classification signal. A linear classifier built solely from these two correlators reproduces approximately 96% of the network predictions (see App. H), even though neither quantity is conventionally regarded as an order parameter of the transition. This illustrates an important distinction between identifying observables conventionally associated with the phase transition and identifying those most useful for a particular learning task. The latter need not coincide with the order parameter or with observables directly tied to the underlying symmetry breaking. For example, Ref. [65] notes that systematic measurement errors lead to a small but nonzero value of ⟨mX ⟩ in the studied dataset. If such a feature were sufficiently predictive of the target labels, TetrisCNN could exploit it. A standard black-box neural network could exploit the same artifact while leaving its role hidden, so high classification accuracy alone could be mistaken for evidence that the network had learned the underlying physics. By contrast, TetrisCNN exposes the observables used by the classifier, allowing them to be scrutinized against theoretical expectations and known experimental imperfections. Physical insight therefore remains necessary to determine whether the identified observables reflect the conventional order parameter or symmetry-breaking physics, other interesting physical structure in the data, or experimental artifacts. d. Task dependence of the learned correlators Moreover, the correlators selected by TetrisCNN depend on

12 the learning objective. In this work, we focused on supervised phase classification, which favors observables that change rapidly across the transition and therefore provide maximal discriminative power between phases. As a consequence, the identified descriptors often resemble conventional phase indicators. However, when the learning objective changes, the network may prefer a different set of correlators better suited to the task at hand. For example, when TetrisCNN is used within prediction-divergence method to regress the experimental tuning parameters from snapshots, it selects correlators informative of continuous variation of the target along the experimental sweep (see App. F 6). e. Rashomon effect Another interesting observation is that the sparsity regularization of TetrisCNN acts primarily as a model-selection mechanism rather than as a performance-enhancing regularizer. As shown in App. F 1, many distinct internal representations achieve nearly identical classification accuracies. Increasing the sparsity regularization progressively removes more complex branches while leaving the predictive performance essentially unchanged, selecting a simple representative from a large family of functionally equivalent solutions. The phenomenon of multiple equally good predictive models existing for the same dataset is known as the Rashomon effect [42] (named after Akira Kurosawa’s movie Rashomon, in which several characters give mutually incompatible accounts of the same event, each plausible from their perspective). It is especially relevant to experimental datasets, as the number of such equally good predictive models tends to increase with data noise. Fortunately for interpretability enthusiasts such as ourselves, the large number of such models also correlates well with the existence of simple yet accurate models [104, 105]. From a physics perspective, it suggests that phase transitions can often be detected in experimental data through various observables. Indeed, several correlators in Fig. 2(b),(d) can be used to distinguish the phases equally well, explaining why different branch combinations can yield nearly identical transition estimates. TetrisCNN, therefore, does not identify a unique set of descriptors relevant to phase transition detection, but rather a particularly simple and interpretable one.

RELATED WORK

Automatic detection of phase transitions and the corresponding order parameters has been a widely recognized goal of machine learning applied to (quantum) physics since 2017 [52]. Rather than replacing existing methods, we view TetrisCNN as complementing the growing toolbox of interpretable approaches. One of the most successful solutions has been to pair neural networks with a renormalization group procedure based on mutual information [106–109], which thrives in the regime of large classical systems and can even be used to study inhomogeneous and biological systems [110].

Tensorial-kernel support vector machines [22, 35, 53–57, 111], thanks to their impressive interpretability, have led to one of the few scientific discoveries with machine learning in physics [22] but at the expense of substantial computational cost and the need to handcraft candidate correlators. Correlator CNN and its extensions [59, 61, 112, 113] learn physically meaningful filters and can even tackle long-range order [61] but require visual inspection to interpret them. Probabilistic variational autoencoders can successfully recover phase structure and candidate phenomena from unlabeled quantum data, including noisy experimental snapshots and nonequilibrium protocols, by learning their compact representation, but interpreting this representation and extracting explicit descriptors still require post-hoc analysis via symbolic regression or substantial physical guidance [33, 34, 84]. Last but not least, one can also resort to exhaustive searches over candidate correlators [114]. TetrisCNN combines several of these desirable properties and offers new ones. It evaluates many candidate correlators simultaneously, automatically selects a sparse subset via regularization, and expresses the resulting classifier directly as symbolic formulas that describe decision boundaries acting on physically meaningful observables. In contrast to post-hoc symbolic-regression approaches [34], whose recovered expressions depend on the userdefined search space and can extrapolate poorly when fitted to the full network output (see App. I), the physical meaning of the TetrisCNN latent coordinates follows directly from the architecture and Boolean-Fourier expansion. The remaining low-dimensional decision rule can then be approximated with a controlled polynomial surrogate. At the same time, the current implementation loses its full interpretability when detecting long-range correlations, allowing only the smallest task-relevant scale to be determined.

VI.

CONCLUSION AND OUTLOOK

Physical models are often sparse descriptions of complex phenomena [115]. In the context of phases of matter, order parameters provide a canonical example of this principle: they provide sparse representations of many-body data sufficient to characterize a phase transition. In this work, we translated this principle into an interpretable convolutional neural-network architecture called TetrisCNN, which uses sparsity as an inductive bias to identify a small set of spin correlators sufficient for solving a given learning task. TetrisCNN combines parallel convolutional branches with filters of different sizes and shapes and promotes sparse bottleneck representations through L1regularization, which favors branches with smaller filters over larger ones. Thanks to the exact mapping between the bottleneck activations and spin correlators, both the latent representation and the learned decision boundaries can be expressed as symbolic formulas in terms of physically meaningful observables. Applied to experimental

13 Rydberg snapshot data for the Ising and XY models, this approach reveals which spin correlators are used to make the prediction and how they are combined to distinguish phases of matter. In the systems studied here, TetrisCNN uncovered crossover and phase indicators that agree with existing theoretical intuition. In the transverse-field Ising dataset, the network bases its prediction on the Z-magnetization, revealing that it detects the crossover out of the trivial polarized state rather than the later emergence of antiferromagnetic order. In the XY dataset, the network identifies a sparse combination of nearest-neighbor and diagonal two-site Z correlators, together with a quadratic contribution from the X-magnetization, closely related to the established in-plane ferromagnetic order parameter. More generally, the ability to express the network predictions in terms of physical observables enables them to be scrutinized, validated, or challenged by physicists. We view TetrisCNN as a step toward machine-assisted discovery. As quantum simulators continue to explore increasingly complex quantum systems, interpretable networks can quickly provide sparse representations of data and help in understanding their new behaviors. To make this happen with TetrisCNN, the next step is to move beyond local order parameters by enriching the branch dictionary to capture nonlocal, string, or other more exotic correlators. Another challenge is to move beyond regular lattices, where geometry itself may generate unfamiliar phases of matter and where standard order-parameter intuition is less developed. More broadly, the TetrisCNN architecture is independent of the learning objective. In this work, we focused on supervised phase classification, but

[1] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017). [2] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Superconducting qubits: Current state of play, Annu. Rev. Condens. Matter Phys. 11, 1 (2019), 1905.13641. [3] A. Browaeys and T. Lahaye, Many-body physics with individually controlled Rydberg atoms, Nat. Phys. 16, 132 (2020), 2002.07413. [4] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021). [5] D. Barredo, V. Lienhard, S. d. Léséleuc, T. Lahaye, and A. Browaeys, Synthetic three-dimensional atomic structures assembled atom by atom, Nature 561, 79 (2018), 1712.02727. [6] J. Carrasquilla, Machine learning for quantum matter, Adv. Phys.: X 5, 1797528 (2020). [7] M. Krenn, J. Landgraf, T. Foesel, and F. Marquardt, Artificial intelligence and machine learning for quantum technologies, Phys. Rev. A 107, 010101 (2023).

the same interpretable architecture can be combined with a wide range of machine-learning paradigms, including distance-learning objectives, prediction-divergence methods, neural quantum states [116], and generative models [84, 117, 118]. This flexibility enables a broad class of machine-learning approaches to become interpretable, extending their role from predictive tools to sources of physical insight into complex many-body systems. Ultimately, we envision machine learning evolving from a tool that makes predictions into one that explains them in the language of physics, enabling an iterative dialogue between experiment, theory, and data-driven models.

DATA AND CODE AVAILABILITY

The data that support the findings of this article and the code used in this study are publicly available at Ref. [66].

ACKNOWLEDGMENTS

We thank Gorka Muñoz-Gil for useful discussions. BvZ and AD acknowledge support from the Dutch National Growth Fund (NGF), as part of the Quantum Delta NL programme and the Top Talent award. The Flatiron Institute is a division of the Simons Foundation. This research was supported in part by grant no. NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). This work was performed using the ALICE compute resources provided by Leiden University.

[8] M. Medvidović and J. R. Moreno, Neural-network quantum states for many-body physics, Eur. Phys. J. Plus 139, 631 (2024). [9] A. Dawid, J. Arnold, B. Requena, A. Gresch, M. Płodzień, K. Donatella, K. A. Nicoli, P. Stornati, R. Koch, M. Büttner, R. Okuła, G. Muñoz-Gil, R. A. Vargas-Hernández, A. Cervera-Lierta, J. Carrasquilla, V. Dunjko, M. Gabrié, P. Huembeli, E. v. Nieuwenburg, F. Vicentini, L. Wang, S. J. Wetzel, G. Carleo, E. Greplová, R. Krems, F. Marquardt, M. Tomza, M. Lewenstein, and A. Dauphin, Machine Learning in Quantum Sciences (Cambridge University Press, 2025). [10] H. Schlömer and A. Bohrdt, Machine learning applications in cold atom quantum simulators (2025), arXiv:2509.08011 [cond-mat.quant-gas]. [11] A. Dawid and Y. LeCun, Introduction to latent variable energy-based models: a path toward autonomous machine intelligence, J. Stat. Mech.: Theory Exp. 2024 (10), 104011, 2306.02572. [12] J. Carrasquilla and R. G. Melko, Machine learning phases of matter, Nat. Phys. 13, 431 (2017). [13] E. P. L. Van Nieuwenburg, Y.-H. Liu, and S. D. Huber, Learning phase transitions by confusion, Nat. Phys. 13, 435 (2017).

14 [14] S. J. Wetzel, Unsupervised learning of phase transitions: From principal component analysis to variational autoencoders, Phys. Rev. E 96, 022140 (2017). [15] E. Greplova, A. Valenti, G. Boschung, F. Schäfer, N. Lörch, and S. D. Huber, Unsupervised identification of topological phase transitions using predictive models, New J. Phys. 22, 045003 (2020). [16] A. Bohrdt, S. Kim, A. Lukin, M. Rispoli, R. Schittko, M. Knap, M. Greiner, and J. Léonard, Analyzing nonequilibrium quantum states through snapshots with artificial neural networks, Phys. Rev. Lett. 127, 150504 (2021). [17] Z. Patel, E. Merali, and S. J. Wetzel, Unsupervised learning of Rydberg atom array phase diagram with siamese neural networks, New J. Phys. 24, 113021 (2022). [18] J. Arnold, F. Schäfer, A. Edelman, and C. Bruder, Mapping out phase diagrams with generative classifiers, Phys. Rev. Lett. 132, 207301 (2024). [19] H. Kim, Y. Zhou, Y. Xu, K. Varma, A. H. Karamlou, I. T. Rosen, J. C. Hoke, C. Wan, J. P. Zhou, W. D. Oliver, Y. D. Lensky, K. Q. Weinberger, and E.-A. Kim, Attention to quantum complexity, Sci. Adv. 11, 41 (2025). [20] O. Malyshev, S. M. Linsel, F. Grusdt, A. Bohrdt, E. Demler, and I. Morera, Distance learning from projective measurements as an information-geometric probe of manybody physics (2026), arXiv:2603.13485 [quant-ph]. [21] K. Kottmann, P. Huembeli, M. Lewenstein, and A. Acín, Unsupervised phase discovery with deep anomaly detection, Phys. Rev. Lett. 125, 170603 (2020). [22] K. Liu, N. Sadoune, N. Rao, J. Greitemann, and L. Pollet, Revealing the phase diagram of Kitaev materials by machine learning: Cooperation and competition between spin liquids, Phys. Rev. Research 3, 023016 (2021). [23] T. Barthel and J. Lu, Fundamental limitations for measurements in quantum many-body systems, Phys. Rev. Lett. 121, 080406 (2018), 1802.04378. [24] R. Ziv, D. Wei, A. Rubio-Abadal, D. Adler, A. Keselman, E. Lustig, R. Talmon, J. Zeiher, I. Bloch, and M. Segev, Unsupervised machine learning for experimental detection of quantum-many-body phase transitions (2025), arXiv:2512.01091 [quant-ph]. [25] B. S. Rem, N. Käming, M. Tarnowski, L. Asteria, N. Fläschner, C. Becker, K. Sengstock, and C. Weitenberg, Identifying quantum phase transitions using artificial neural networks on experimental data, Nat. Phys. 15, 917 (2019). [26] E. Khatami, E. Guardado-Sanchez, B. M. Spar, J. F. Carrasquilla, W. S. Bakr, and R. T. Scalettar, Visualizing strange metallic correlations in the two-dimensional fermi-hubbard model with artificial intelligence, Phys. Rev. A 102, 033326 (2020). [27] N. Käming, A. Dawid, K. Kottmann, M. Lewenstein, K. Sengstock, A. Dauphin, and C. Weitenberg, Unsupervised machine learning of topological phase transitions from experimental data, Mach. Learn.: Sci. Technol. 2, 035037 (2021). [28] 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 (2023). [29] C. Miles, R. Samajdar, S. Ebadi, T. T. Wang, H. Pichler, S. Sachdev, M. D. Lukin, M. Greiner, K. Q. Weinberger, and E.-A. Kim, Machine learning discovery of new phases in programmable quantum simulator snapshots, Phys.

Rev. Res. 5, 013026 (2023). [30] H. Zhang, S. Jiang, X. Wang, W. Zhang, X. Huang, X. Ouyang, Y. Yu, Y. Liu, D.-L. Deng, and L.-M. Duan, Experimental demonstration of adversarial examples in learning topological phases, Nat. Commun. 13, 4993 (2022), 2111.12715. [31] A. Madry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu, Towards deep learning models resistant to adversarial attacks, in International Conference on Learning Representations (2018). [32] Y. Yu, L.-W. Yu, W. Zhang, H. Zhang, X. Ouyang, Y. Liu, D.-L. Deng, and L.-M. Duan, Experimental unsupervised learning of non-hermitian knotted phases with solid-state spins, Npj Quantum Inf. 8, 116 (2022). [33] P. d. Schoulepnikoff, G. Muñoz-Gil, H. P. Nautrup, and H. J. Briegel, Interpretable representation learning of quantum data enabled by probabilistic variational autoencoders, Phys. Rev. A 112, 062423 (2025). [34] P. d. Schoulepnikoff, H. P. Nautrup, H. J. Briegel, and G. Muñoz-Gil, Discovering quantum phenomena with interpretable machine learning (2026), arXiv:2604.16015 [quant-ph]. [35] N. Sadoune, I. Pogorelov, C. L. Edmunds, G. Giudici, G. Giudice, C. D. Marciniak, M. Ringbauer, T. Monz, and L. Pollet, Learning symmetry-protected topological order from trapped-ion experiments, Quantum 10, 2100 (2026), 2408.05017. [36] N. Elhage, T. Hume, C. Olsson, N. Schiefer, T. Henighan, S. Kravec, Z. Hatfield-Dodds, R. Lasenby, D. Drain, C. Chen, R. Grosse, S. McCandlish, J. Kaplan, D. Amodei, M. Wattenberg, and C. Olah, Toy models of superposition (2022). [37] D. Klindt, C. O’Neill, P. Reizinger, H. Maurer, and N. Miolane, From superposition to sparse codes: interpretable representations in neural networks (2025), arXiv:2503.01824 [cs.LG]. [38] K. Vafa, J. Y. Chen, A. Rambachan, J. Kleinberg, and S. Mullainathan, Evaluating the world model implicit in a generative model, in Advances in Neural Information Processing Systems, Vol. 37 (Curran Associates, Inc., 2024) pp. 26941–26975. [39] K. Vafa, P. G. Chang, A. Rambachan, and S. Mullainathan, What has a foundation model found? using inductive bias to probe for world models, in Proceedings of the 42nd International Conference on Machine Learning, ICML’25, Vol. 267 (2025) pp. 60727–60747. [40] A. Ghorbani, A. Abid, and J. Zou, Interpretation of neural networks is fragile, Proceedings of the AAAI Conference on Artificial Intelligence 33, 3681 (2019). [41] K. Cybiński, M. Płodzień, M. Tomza, M. Lewenstein, A. Dauphin, and A. Dawid, Characterizing out-ofdistribution generalization of neural networks: application to the disordered Su–Schrieffer–Heeger model, Mach. Learn.: Sci. Technol. 6, 015014 (2025). [42] C. Rudin, C. Zhong, L. Semenova, M. Seltzer, R. Parr, J. Liu, S. Katta, J. Donnelly, H. Chen, and Z. Boner, Position: Amazing things come from having many good models, in Proceedings of the 41st International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 235, edited by R. Salakhutdinov, Z. Kolter, K. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp (PMLR, 2024) pp. 42783–42795. [43] A. Dawid, P. Huembeli, M. Tomza, M. Lewenstein, and A. Dauphin, Phase detection with neural networks: in-

15 terpreting the black box, New J. Phys. 22, 115001 (2020). [44] S. J. Wetzel, R. G. Melko, J. Scott, M. Panju, and V. Ganesh, Discovering symmetry invariants and conserved quantities by interpreting siamese neural networks, Phys. Rev. Res. 2, 033499 (2020). [45] A. Dawid, P. Huembeli, M. Tomza, M. Lewenstein, and A. Dauphin, Hessian-based toolbox for reliable and interpretable machine learning in physics, Mach. Learn.: Sci. Techn. 3, 015002 (2021). [46] J. Arnold, F. Schäfer, M. Žonda, and A. U. J. Lode, Interpretable and unsupervised phase classification, Phys. Rev. Research 3, 033052 (2021). [47] J. Arnold and F. Schäfer, Replacing neural networks by optimal analytical predictors for the detection of phase transitions, Phys. Rev. X 12, 031044 (2022). [48] J. Arnold, N. Lörch, F. Holtorf, and F. Schäfer, Machine learning phase transitions: Connections to the Fisher information (2023), arXiv:2311.10710 [cond-mat.dis-nn]. [49] S. J. Wetzel, Closed-form interpretation of neural network classifiers with symbolic regression gradients (2024), arXiv:2401.04978 [cs.LG]. [50] Y. Zhan, A. Elben, H.-Y. Huang, and Y. Tong, Learning conservation laws in unknown quantum dynamics, PRX Quantum 5, 010350 (2024). [51] Y. Ju, S. S. Alam, J. Minoff, F. Anselmi, H. Pu, and A. Patel, Interpreting convolutional neural networks’ low-dimensional approximation to quantum spin systems, Phys. Rev. Research 7, 013094 (2025). [52] S. J. Wetzel and M. Scherzer, Machine learning of explicit order parameters: From the Ising model to SU(2) lattice gauge theory, Phys. Rev. B 96, 184410 (2017). [53] J. Greitemann, K. Liu, and L. Pollet, Probing hidden spin order with interpretable machine learning, Phys. Rev. B 99, 060404 (2019). [54] K. Liu, J. Greitemann, and L. Pollet, Learning multiple order parameters with interpretable machines, Phys. Rev. B 99, 104410 (2019). [55] J. Greitemann, K. Liu, L. D. C. Jaubert, H. Yan, N. Shannon, and L. Pollet, Identification of emergent constraints and hidden order in frustrated magnets using tensorial kernel methods of machine learning, Phys. Rev. B 100, 174408 (2019). [56] N. Sadoune, G. Giudici, K. Liu, and L. Pollet, Unsupervised interpretable learning of phases from many-qubit systems, Phys. Rev. Res. 5, 013082 (2023). [57] N. Sadoune, K. Liu, H. Yan, L. D. C. Jaubert, N. Shannon, and L. Pollet, Human-machine collaboration: Ordering mechanism of rank-2 spin liquid on breathing pyrochlore lattice, Phys. Rev. Research 7, 033061 (2025), 2402.10658. [58] A. Cole, G. J. Loges, and G. Shiu, Quantitative and interpretable order parameters for phase transitions from persistent homology, Phys. Rev. B 104, 104426 (2021). [59] C. Miles, A. Bohrdt, R. Wu, C. Chiu, M. Xu, G. Ji, M. Greiner, K. Q. Weinberger, E. Demler, and E.-A. Kim, Correlator convolutional neural networks as an interpretable architecture for image-like quantum matter data, Nat. Commun. 12, 3905 (2021). [60] S. Striegel, E. Ibarra-García-Padilla, and E. Khatami, Machine learning detection of correlations in snapshots of ultracold atoms in optical lattices (2023), arXiv:2310.03267 [cond-mat.dis-nn]. [61] H. Schlömer and A. Bohrdt, Fluctuation based interpretable analysis scheme for quantum many-body snap-

shots, SciPost Phys. 15, 099 (2023). [62] C. Cao, F. M. Gambetta, A. Montanaro, and R. A. Santos, Unveiling quantum phase transitions from traps in variational quantum algorithms, Npj Quantum Inf. 11, 93 (2025). [63] A. Suresh, H. Schlömer, B. Hashemi, and A. Bohrdt, Interpretable correlator transformer for image-like quantum matter data, Mach. Learn.: Sci. Technol. 6, 025006 (2025). [64] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum simulation of 2D antiferromagnets with hundreds of Rydberg atoms, Nature 595, 233–238 (2021). [65] C. Chen, G. Bornet, M. Bintz, G. Emperauger, L. Leclerc, V. S. Liu, P. Scholl, D. Barredo, J. Hauschild, S. Chatterjee, M. Schuler, A. M. Läuchli, M. P. Zaletel, T. Lahaye, N. Y. Yao, and A. Browaeys, Continuous symmetry breaking in a two-dimensional Rydberg array, Nature 616, 691–695 (2023). [66] K. Cybiński, B. van Zwol, J. Enouen, G. Bornet, T. Lahaye, A. Browaeys, A. Georges, and A. Dawid, https: //doi.org/10.5281/zenodo.14035852 (2026), GitHub repository: TetrisCNN for spin systems (Version arXiv 2.0). [67] H. Labuhn, D. Barredo, S. Ravets, S. d. Léséleuc, T. Macrì, T. Lahaye, and A. Browaeys, Tunable twodimensional arrays of single Rydberg atoms for realizing quantum ising models, Nature 534, 667 (2016), 1509.04543. [68] O. Simard, A. Dawid, J. Tindall, M. Ferrero, A. M. Sengupta, and A. Georges, Learning interactions between Rydberg atoms, PRX Quantum 6, 030324 (2025), 2412.12019. [69] K. Cybinski, J. Enouen, A. Georges, and A. Dawid, Speak so a physicist can understand you! TetrisCNN for detecting phase transitions and order parameters, in Machine Learning and the Physical Sciences Workshop @ NeurIPS 2024 (2024). [70] L. D. Landau, On the theory of phase transitions, Zh. Eksp. Teor. Fiz. 7, 19 (1937). [71] P. Hohenberg and A. Krekhov, An introduction to the ginzburg–landau theory of phase transitions and nonequilibrium patterns, Phys. Rep. 572, 1 (2015). [72] V. L. Berezinskii, Destruction of long-range order in one-dimensional and two-dimensional systems having a continuous symmetry group i. classical systems, Sov. Phys. JETP 32, 493 (1971). [73] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, J. Phys. C: Solid State Phys. 6, 1181 (1973). [74] X.-G. Wen, Topological orders in rigid states, Int. J. Mod. Phys. B 4, 239 (1990). [75] X.-G. Wen and Q. Niu, Ground-state degeneracy of the fractional quantum Hall states in the presence of a random potential and on high-genus Riemann surfaces, Phys. Rev. B 41, 9377 (1990). [76] A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Ann. Phys. 303, 2 (2003), quant-ph/9707021. [77] M. d. Nijs and K. Rommelse, Preroughening transitions in crystal surfaces and valence-bond phases in quantum spin chains, Phys. Rev. B 40, 4709 (1989). [78] T. Hikihara, L. Kecke, T. Momoi, and A. Furusaki, Vector chiral and multipolar orders in the spin- 12 frustrated

16 ferromagnetic chain in magnetic field, Phys. Rev. B 78, 144404 (2008), 0807.0858. [79] M. Kalinowski, R. Samajdar, R. G. Melko, M. D. Lukin, S. Sachdev, and S. Choi, Bulk and boundary quantum phase transitions in a square Rydberg atom array, Phys. Rev. B 105, 174417 (2022), 2112.10790. [80] M. J. O’Rourke and G. K.-L. Chan, Entanglement in the quantum phases of an unfrustrated Rydberg atom array, Nat. Commun. 14, 5397 (2023), 2201.03189. [81] L. Wang, Discovering phase transitions with unsupervised learning, Phys. Rev. B 94, 195105 (2016). [82] A. Lidiak and Z. Gong, Unsupervised machine learning of quantum phase transitions using diffusion maps, Phys. Rev. Lett. 125, 225701 (2020), 2003.07399. [83] K. Kottmann, F. Metz, J. Fraxanet, and N. Baldelli, Variational quantum anomaly detection: Unsupervised mapping of phase diagrams on a physical quantum computer, Phys. Rev. Research 3, 043184 (2021), 2106.07912. [84] F. Møller, G. Fernández-Fernández, T. Schweigler, P. d. Schoulepnikoff, J. Schmiedmayer, and G. Muñoz-Gil, Learning minimal representations of many-body physics from snapshots of a quantum simulator, Phys. Rev. Research 8, 023094 (2026), 2509.13821. [85] F. Schäfer and N. Lörch, Vector field divergence of predictive model output as indication of phase transitions, Phys. Rev. E 99, 062107 (2019). [86] T. Chalopin, I. Ferrier-Barbut, T. Lahaye, A. Browaeys, and D. Clément, Connected correlations in cold atom experiments, Comptes Rendus. Physique 27, 65 (2026). [87] R. Eldan and O. Shamir, The power of depth for feedforward neural networks, in 29th Annual Conference on Learning Theory, Vol. 49, edited by V. Feldman, A. Rakhlin, and O. Shamir (PMLR, Columbia University, New York, New York, USA, 2016) pp. 907–940. [88] J. Kaplan, S. McCandlish, T. Henighan, T. B. Brown, B. Chess, R. Child, S. Gray, A. Radford, J. Wu, and D. Amodei, Scaling laws for neural language models (2020), arXiv:2001.08361 [cs.LG]. [89] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, Nature 521, 436 (2015). [90] A. G. Wilson, Position: Deep learning is not so mysterious or different, in Proceedings of the 42nd International Conference on Machine Learning, Vol. 267 (2025) pp. 82326–82346. [91] R. Agarwal, L. Melnick, N. Frosst, X. Zhang, B. Lengerich, R. Caruana, and G. E. Hinton, Neural additive models: Interpretable machine learning with neural nets, in Advances in Neural Information Processing Systems, Vol. 34 (Curran Associates, Inc., 2021) pp. 4699–4711. [92] J. Enouen and Y. Liu, Sparse interaction additive networks via feature interaction detection and sparse selection, in Proceedings of the 36th International Conference on Neural Information Processing Systems, NIPS ’22 No. 1011 (Curran Associates Inc., Red Hook, NY, USA, 2022) pp. 13908–1392. [93] C. Rudin, C. Chen, Z. Chen, H. Huang, L. Semenova, and C. Zhong, Interpretable machine learning: Fundamental principles and 10 grand challenges, Statistic Surveys 16, 1 (2022). [94] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, Gradient-based learning applied to document recognition, in Proceedings of the IEEE , Vol. 86 (1998) pp. 2278– 2324.

[95] T. Cohen and M. Welling, Group equivariant convolutional networks, in Proceedings of The 33rd International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 48 (PMLR, 2016) pp. 2990–2999. [96] R. O’Donnell, Analysis of Boolean Functions (Cambridge University Press, 2014). [97] I. Schurov, A. Kravchenko, M. I. Katsnelson, A. A. Bagrov, and T. Westerhout, Learning complexity of many-body quantum sign structures through the lens of Boolean Fourier analysis (2025), arXiv:2508.09870 [cond-mat.dis-nn]. [98] F. Döschl and A. Bohrdt, Towards interpretability of neural quantum states (2025), arXiv:2508.14152 [quantph]. [99] M. Nicolau, A. R. Tavares, Z. Zhang, P. H. C. Avelar, J. M. Flach, L. D. Lamb, and M. Vardi, Understanding Boolean function learnability on deep neural networks: PAC learning meets neurosymbolic models, in Proceedings of The 19th International Conference on Neurosymbolic Learning and Reasoning, Vol. 284 (PMLR, 2025) pp. 719–735. [100] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner, betaVAE: Learning basic visual concepts with a constrained variational framework, in ICLR 2017 - 5th International Conference on Learning Representations (2017) arXiv:1804.03599. [101] H. Kim and A. Mnih, Disentangling by factorising, in Proceedings of the 35th International Conference on Machine Learning, Vol. 80 (PMLR, 2018) pp. 2649–2658. [102] H. Cunningham, A. Ewart, L. Riggs, R. Huben, and L. Sharkey, Sparse autoencoders find highly interpretable features in language models, in ICLR 2024 - 12th International Conference on Learning Representations (2024). [103] T. Mendes-Santos, M. Schmitt, A. Angelone, A. Rodriguez, P. Scholl, H. J. Williams, D. Barredo, T. Lahaye, A. Browaeys, M. Heyl, and M. Dalmonte, Wave-function network description and Kolmogorov complexity of quantum many-body systems, Phys. Rev. X 14, 021029 (2024). [104] L. Semenova, C. Rudin, and R. Parr, On the existence of simpler machine learning models, 2022 ACM Conference on Fairness Accountability and Transparency , 1827 (2022), 1908.01755. [105] L. Semenova, H. Chen, R. Parr, and C. Rudin, A path to simpler models starts with noise, in Thirty-seventh Conference on Neural Information Processing Systems (2023). [106] M. Koch-Janusz and Z. Ringel, Mutual information, neural networks and the renormalization group, Nat. Phys. 14, 578 (2018), 1704.06279. [107] D. E. Gökmen, Z. Ringel, S. D. Huber, and M. KochJanusz, Symmetries and phase diagrams with real-space mutual information neural estimation, Phys. Rev. E 104, 064106 (2021), 2103.16887. [108] D. E. Gökmen, Z. Ringel, S. D. Huber, and M. KochJanusz, Statistical physics through the lens of real-space mutual information, Phys. Rev. Lett. 127, 240603 (2021), 2101.11633. [109] A. Gordon, A. Banerjee, M. Koch-Janusz, and Z. Ringel, Relevance in the renormalization group and in information theory, Phys. Rev. Lett. 126, 240601 (2021), 2012.01447.

17 [110] D. E. Gökmen, S. Biswas, S. D. Huber, Z. Ringel, F. Flicker, and M. Koch-Janusz, Compression theory for inhomogeneous systems, Nat. Commun. 15, 10214 (2024), 2301.11934. [111] N. Rao, K. Liu, and L. Pollet, Inferring hidden symmetries of exotic magnets from detecting explicit order parameters, Phys. Rev. E 104, 015311 (2021), 2007.07000. [112] Z. Gong, A. Halder, A. Bohrdt, S. Seitz, and D. Gebauer, C3nn: Cosmological correlator convolutional neural network an interpretable machine-learning framework for cosmological analyses, Astrophys. J. 971, 156 (2024). [113] A. Suresh, H. Schlömer, B. Hashemi, and A. Bohrdt, Interpretable correlator transformer for image-like quantum matter data, Mach. Learn.: Sci. Technol. 6, 025006 (2025). [114] R. Verdel, V. Vitale, R. K. Panda, E. D. Donkor, A. Rodriguez, S. Lannig, Y. Deller, H. Strobel, M. K. Oberthaler, and M. Dalmonte, Data-driven discovery of statistically relevant information in quantum simulators, Phys. Rev. B 109, 075152 (2024). [115] D. A. Roberts, Why is AI hard and physics simple? (2021), arXiv:2104.00008 [hep-th]. [116] A. Valenti, E. Greplova, N. H. Lindner, and S. D. Huber, Correlation-enhanced neural networks as interpretable variational quantum states, Phys. Rev. Research 4, L012010 (2022), 2103.05017. [117] C. Casert, K. Mills, T. Vieijra, J. Ryckebusch, and I. Tamblyn, Optical lattice experiments at unobserved conditions with generative adversarial deep learning, Phys. Rev. Research 3, 033267 (2021). [118] D. Fitzek, Y. H. Teoh, H. P. C. Fung, G. A. Dagnew, E. Merali, M. S. Moss, B. MacLellan, and R. G. Melko, RydbergGPT, Mach. Learn.: Sci. Technol. 6, 045057 (2025). [119] J. Gross, Combinatorial Methods with Computer Applications, Discrete Mathematics and Its Applications (Taylor & Francis, 2008). [120] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, in International Conference on Learning Representations (2015). [121] I. Loshchilov and F. Hutter, Decoupled weight decay regularization, in International Conference on Learning Representations (2019). [122] L. McInnes, J. Healy, and J. Melville, UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction (2020), arXiv:1802.03426 [stat.ML]. [123] R. J. G. B. Campello, D. Moulavi, and J. Sander, Density-based clustering based on hierarchical density estimates, in Advances in Knowledge Discovery and Data Mining (Springer Berlin Heidelberg, 2013) p. 160–172. [124] J. Ansel, E. Yang, H. He, N. Gimelshein, A. Jain, M. Voznesensky, B. Bao, P. Bell, D. Berard, E. Burovski, G. Chauhan, A. Chourdia, W. Constable, A. Desmaison, Z. DeVito, E. Ellison, W. Feng, J. Gong, M. Gschwind, B. Hirsh, S. Huang, K. Kalambarkar, L. Kirsch, M. Lazos, M. Lezcano, Y. Liang, J. Liang, Y. Lu, C. K. Luk, B. Maher, Y. Pan, C. Puhrsch, M. Reso, M. Saroufim, M. Y. Siraichi, H. Suk, S. Zhang, M. Suo, P. Tillet, X. Zhao, E. Wang, K. Zhou, R. Zou, X. Wang, A. Mathews, W. Wen, G. Chanan, P. Wu, and S. Chintala, Pytorch 2: Faster machine learning through dynamic python bytecode transformation and graph compilation, in Proceedings of the 29th ACM International Conference on Architectural Support for Programming Languages and

Operating Systems, Volume 2 , ASPLOS ’24 (Association for Computing Machinery, New York, NY, USA, 2024) p. 929–947. [125] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Scikitlearn: Machine learning in Python, J. Mach. Learn. Res. 12, 2825 (2011). [126] M. T. Ribeiro, S. Singh, and C. Guestrin, “Why Should I Trust You?”: Explaining the predictions of any classifier, in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (2016) pp. 1135–1144. [127] M. Cranmer, Interpretable machine learning for science with PySR and SymbolicRegression.jl (2023), arXiv:2305.01582 [astro-ph.IM].

Contents of the appendices A. Data preparation A 1. Snapshot imbalance and data splitting A 2. Combining two measurement bases B. Correlators B 1. Boolean function theory B 2. Counting patterns C. Details of TetrisCNN architecture and hyperparameters D. Separation of active and deactivated branches by learning-rate annealing D 1. Origin of the activation floor D 2. Active versus deactivated branches D 3. Sharpening sparsity with learning-rate annealing E. Unsupervised detection of transitions in the experimental sweeps: method comparison E 1. Learning by confusion E 2. Prediction-divergence method E 3. Dimensionality reduction and clustering E 4. Distance learning E 5. Comparison of transition estimates F. Supporting material F 1. Robustness of the learned descriptors to regularization F 2. Searching for the relevant scale F 3. Rotationally invariant and unconstrained TetrisCNN F 4. Robustness to random X/Z pairings F 5. Contribution of the X and Z bases to classification F 6. Task dependence of the learned descriptors F 7. Complete mapping from filters to correlators G. R-squared for multivariate and aggregate output H. Fitting the decision boundary to the classifier predictions in the bottleneck space H 1. Fitting the decision boundary H 2. Decision boundary: plane vs. quadratic surface I. Limitations of symbolic regression

18 Appendix A: Data preparation Snapshot imbalance across the sweep and data splitting

A characteristic challenge of working with raw projective measurements from quantum simulators, as opposed to numerically simulated ground states, is the highly nonuniform distribution of snapshots acquired across the experimental sweep. As shown in Fig. 7, the number of snapshots varies by up to an order of magnitude along the experimental sweep. Throughout this work, we therefore use a stratified 70/30 train-validation split performed independently at each sampling point, ensuring that all regions of the sweep are represented in both sets. We verified that our results are insensitive to this choice by comparing the stratified split with three alternatives: a global random split, a stratified split with the number of snapshots per sampling point capped at 500, and training (a) Number of Snapshots

1500

Total: 20895 snapshots

Min: 297

1250 Across 28 time points

0.0

0.2

0.4

0.6

0.8

12.5

δXY /2π (MHz)

15.0

10.0

1.0 Z δXY X δXY XZ δXY

9 7.5

5.0

2.5 0.0 0

t (µs)

FIG. 8: The two measurement bases are acquired in separate sweeps whose control fields do not coincide. Realized staggered detuning δXY across the experimental sweeps for the independent X-basis and Z-basis measurements of the XY system. The inset highlights the discrepancy during the initial rapid XZ (green ramp-down. The effective tuning parameter δXY curve) is the arithmetic mean of the two physical sweeps, used when combining snapshots from both measurement bases.

Min shared snapshots: 297

1000

with an inverse-frequency-weighted loss. Across these choices, we observe only modest differences in validation accuracy, within 2 percentage points. We therefore retain the stratified split as the simplest approach that preserves all available experimental snapshots without modifying their contribution to the loss.

750 500 250

0. 0.6 1.8 1.0 1.2 1.4 1.6 2.8 2.0 2.2 2.4 2.6 3.8 3.0 3.2 3.4 3.6 4.8 4.0 4.2 4.4 4.6 5.8 5.0 5.2 5.4 5.6 6.8 0

t (µs)

Combining two measurement bases

(b) Number of Snapshots

Total snapshots count:

6000 Z: 26924

X: 16484 XZ: 14471 4000 Across 21 time points Min shared snapshots: 371 5000

Min: 371 Z basis X basis XZ basis (X+Z)

3000 2000 1000

0. 0.0 12 0. 5 0. 25 37 5 0. 0.5 67 0. 5 75 1. 1.0 12 1 5 1. .25 37 5 1. 1.5 67 1 5 1. .75 87 5 2. 0 2. 5 3. 0 4. 0 6. 0 8. 0

t (µs)

FIG. 7: Snapshot counts vary by up to an order of magnitude across the sweep. Distribution of the acquired snapshots over the parameter sweep steps for (a) the 8 × 8 transverse-field Ising model and (b) the 6 × 7 XY model. Red dashed lines indicate the minimum snapshot count (297 for the Ising dataset and 371 for the XY dataset).

Projective measurements are destructive, so the Xand Z-basis snapshots of the XY dataset originate from separate experimental realizations performed under nominally identical conditions. Combining them therefore requires addressing two differences between the datasets: small discrepancies in the realized control parameters and the absence of a ground-truth pairing between individual snapshots. a. Differing Hamiltonian parameters. As shown in Fig. 8, the staggered detuning δXY of the XY experiment is similar between the two measurement sweeps except during the initial rapid ramp-down, where the discrepancy reaches approximately ∼ 3 MHz. For the combined XZ dataset, we therefore define the effective tuning parameter XZ X Z as the arithmetic mean, δXY = δXY + δXY /2. b. Pairing snapshots between bases. Because the Xand Z-basis snapshots are independent experimental realizations, there is no unique pairing between individual snapshots acquired at the same sampling point. We combine the two datasets by pairing snapshots within each sampling point according to their acquisition index, until the smaller dataset is exhausted (Fig. 7). This pairing

19 is therefore arbitrary. In App. F 4, we verify that the selected correlators and recovered decision boundary are robust to random re-pairings of the X- and Z-basis snapshots.

to define a pattern. The number of patterns of all orders on D sites is thus: D X D−1 k=1

k−1

= 2D−1

Appendix B: Correlators 1.

Boolean function theory

The fact that snapshot data takes the form S ∈ {±1}|I| where I is an index set, means that neural networks trained on this data are Boolean functions. A basic result in Boolean function analysis is the Fourier expansion theorem, which states that any f : {±1}|I| → R can be uniquely expressed as a multilinear polynomial: f (S) =

cA (f )

A⊆I

(B1)

j∈A

where the sum is over all subsets of I, including the empty subset. cA are (real-valued) Fourier coefficients Q in the expansion, and j∈A Sj are the basis functions (known in the literature as monomials, characters, or parity functions). As stated in the main text, we focus on convolutions of spatially arranged Boolean data. Following Eq. (5), we define f [P ]T = f [P ](S)T as a function depending on spins {Si }i∈T (P ) within a translated pattern P . With Eq. (B1) (I = P , A = P ′ ), this gives: f [P ]T =

X P ′ ⊆P

cP ′ (f )

i∈T (P ′ )

Applying a spatial average to this gives: Y 1 X 1 X X f [P ]T = cP ′ (f ) Si |TP | |TP | T T P ′ ⊆P i∈T (P ′ ) X = cP ′ (f ) C[P ′ ], P ′ ⊆P

which is the result used in the main text. Note that this requires TP ′ = TP for all P ′ , i.e. correlators C[P ′ ] need to be calculated over the same spins as C[P ] (e.g. C[□■] is only equivalent to C[■] up to edge effects).

Counting patterns

The number of patterns in 1D can be worked out as follows. It amounts to choosing k positions from {1, ..., D} with translation invariance; shifted patterns are equivalent. That is, instead of considering k absolute sites, one can consider instead k − 1 ‘spacings’ (relative distances) between sites. Having fixed the first site (being equivalent to all others), choose k − 1 from the D − 1 remaining sites

For 2D, finding the number of patterns is more complicated due to additional symmetries (rotations and reflections), but the number similarly scales exponentially and can be found using Burnside-Pólya Counting [119]. The numbers in the main text were found instead by complete enumeration with a simple Python script (see Tab. II). TABLE II: Number of subpatterns for 3x3 and 4x4 filters and different symmetry groups. symmetries

3x3

4x4

– 511 65,535 translations 400 57,856 rotations, reflections 101 8,547 translations, rotations, reflections 85 7,625

Appendix C: Details of TetrisCNN architecture and hyperparameters

This Appendix complements the architectural description of TetrisCNN in Sec. II C. A detailed schematic of the architecture is shown in Fig. 9, while the hyperparameters used for the classifiers in the main text are summarized in Tab. III. a. Task network The task neural network is a flexible readout of the interpretable bottleneck and is not intrinsic to TetrisCNN. For the classification and regression tasks considered here, we use a multilayer perceptron. We tested architectures with one to three hidden layers and found no dependence of performance on the task network architecture. We therefore use a twohidden-layer fully-connected network with architecture [ #branches, 32, 16, 2 ] for classification. b. Branch penalty As defined in Eq. (11), the branch penalty λk increases with the pattern cardinality |Pk | through a log-linear interpolation between λmin and λmax , favoring lower-order correlators. More generally, the branch penalty provides a way to encode what is considered a simple physical description. For example, when including dilated patterns, the penalty can increase with spatial extent to favor local correlators. Similarly, nonequivariant branches can be penalized more strongly than their symmetry-constrained counterparts to favor symmetric descriptions. Thus, the penalty can be adapted to encode prior knowledge about the complexity of candidate physical observables.

FIG. 9: Detailed architecture of TetrisCNN. (a) A batch of B input snapshots (or pairs of snapshots, when two measurement bases are combined) is processed by (b) K parallel convolutional branches defined by patterns Pk . Each branch produces one scalar activation zk through two convolutional layers and global average pooling [Eq. (8)]. The first convolution sets the spatial pattern Pk and is the only stage that reads correlations between sites; the 1 × 1 convolution acts site-by-site across channels and therefore cannot introduce any new spatialP correlations. The resulting activations {zk }K form the sparse interpretable bottleneck, on which the L1 penalty k=1 k λk |zk | acts. (c) The bottleneck is then read out by a task network (here, a small fully connected network) whose architecture is adapted to the learning task.

Appendix D: Separation of active and deactivated branches by learning-rate annealing 1.

Origin of the activation floor

Everything we read off the TetrisCNN bottleneck relies on distinguishing active branches from deactivated ones. In practice, deactivated branches do not decay exactly to zero but settle onto a small residual floor, visible as the nearly flat lines at the bottom of Fig. 4(b,e). Here we explain the origin of this floor and why learning-rate annealing is required to clearly separate active and deactivated branches, yielding a sparse bottleneck. Let z ∈ RK denote the bottleneck activations of the K parallel branches for a given input [cf. Eq. (8)]. The network is trained with the task loss together with the branch-wise sparsity penalty, L = LTask (ŷ, y) +

K X

updates induce corresponding changes in zk through the branch Jacobian. We return to this point later. Therefore, differentiating with respect to a single bottleneck activation gives ∇zk L = ∇zk LTask + λk sgn(zk ) . | {z } | {z } task signal

The second term always pushes the activation toward zero. Once the typical task contribution of a branch becomes negligible compared with its sparsity penalty, AvgB [|∇zk LTask |] ≪ λk ,

(D3)

the gradient entering the optimizer is dominated by the L1 term. For an optimization step t using minibatch Bt , we can therefore write (t)

λk |zk | ,

(D1)

k=1

where λk > 0 is the penalty on branch k. To expose the mechanism behind the activation floor, we first treat the bottleneck activations zk as effective optimization variables. In the actual network, AdamW updates the branch parameters rather than zk directly, but parameter

(D2)

sparsity term

gt = ∇zk L(t) ≈ λk sgn(zk ).

(D4)

For ordinary gradient descent, this would give (t)

(t)

∆zk = −αt gt ≈ −αt λk sgn(zk ),

(D5)

so the magnitude of the update would scale directly with λk and instantaneous learning rate αt . Adam or AdamW (which differs from Adam only in how weight decay is

21 TABLE III: Hyperparameters of the two classification models reported in the main text. The columns correspond to the Ising and XY classifiers shown in panels (b) and (f) of Fig. 4, respectively. Entries spanning both columns are shared between the two models. TFIM (Z basis) XY (Z, X bases)

Dataset TetrisCNN

[ ■, stride=1, dilation=1 ] [ ■■, stride=1, dilation=1, rot=C4 ] [ ■□ □■, stride=1, dilation=1, rot=C4 ] [ ■■ ■□, stride=1, dilation=1, rot=C4 ]

Available branches

[ ■■ ■■, stride=1, dilation=1, rot=C4 ] Bottleneck size (# branches) λmin , λmax

5 10−3 , 103

Task NN Task type Architecture Layers

Classification Multilayer Perceptron [# branches, 32, 16, 2]

Training parameters Loss function Optimizer Weight decay Learning rate scheduler Initial learning rate Reduce factor Patience Minimal learning rate Max epochs Early Stopping Warmup Patience RNG seed Batch size

Cross Entropy Loss AdamW 1e-05 Reduce on Plateau 1e-2 0.5 5 1e-9 250 Yes 150 10 42 64

|gt | ≈ λk ,

(D7)

Substituting this into Eq. (D6) gives (t)

λk sgn(zk ) (t) = −αt sgn(zk ). λk

In the full network, AdamW updates the branch parameters θ k rather than the activation zk directly. In the penalty-dominated regime, the gradient of the sparsity term with respect to a branch parameter θk,j is, by the chain rule,

∆zk ≈ ∇θk zk · ∆θ k = O(αt ).

If this penalty-dominated regime persists sufficiently long, Adam’s running second-moment estimate approaches the squared gradient scale, √ vt ∼ λ2k , vt ∼ λk . (D8)

∆zk ∼ −αt

(D10)

(D11)

The corresponding Adam second-moment estimate therefore scales as λ2k (∂zk /∂θk,j )2 . As in the simplified activation-level argument, the overall factor λk cancels between the gradient and its normalization. The resulting parameter update is therefore of order αt , with a prefactor determined by the local derivatives of zk with respect to the branch parameters. For a sufficiently small parameter update,

where vt is the exponential moving average of gt2 . For the penalty-dominated gradient of Eq. (D4), its magnitude is approximately constant, gt2 ≈ λ2k .

zfloor = O(αt ).

∂zk ∂LL1 = λk sgn(zk ) . ∂θk,j ∂θk,j

applied) [120, 121] behave differently because they normalize the gradient by the square root of its running second moment. Ignoring, for the moment, the first-moment dynamics, the Adam update has the schematic form gt (t) ∆zk ∼ −αt √ , (D6) vt

(t)

whether the branch enters this penalty-dominated regime, but once it does, increasing λk no longer substantially reduces the residual update scale under Adam. Near zero, these updates repeatedly change the sign of zk . For example, if 0 < zk < αt , the next update can carry it across zero. The sign of the L1 gradient then reverses and the following updates push it back in the opposite direction. The activation therefore fluctuates around zero instead of converging smoothly to it. The characteristic scale of these residual fluctuations is set by the instantaneous learning rate,

(D9)

The important point is that the leading dependence on λk cancels. The sparsity strength still determines

(D12)

Thus, moving from the simplified activation-level picture to the actual parameter updates changes the precise prefactor but preserves the central result: once the sparsity term dominates, the residual activation scale is controlled by the learning rate rather than the overall magnitude of λk .

Active versus deactivated branches

A branch remains active when its contribution to the task is sufficiently large to counteract the sparsity penalty. Near convergence, this requires the task and sparsity contributions to approximately balance, AvgB [|∇zk LTask |] ∼ λk .

(D13)

If instead the task signal becomes much smaller than λk , the branch is driven toward the learning-rate-dependent floor described in Sec. D 1. This is the mechanism by which the branch penalties perform model selection. Increasing λmax raises the penalty on the more complex patterns, so branches whose contribution to the task is insufficient are progressively suppressed. In practice, the surviving set rapidly becomes

22 stable over a broad range of λmax while the classification accuracy remains essentially unchanged, as shown in Fig. 15 of App. F 1. We therefore identify a branch as active when its converged activation remains well above the residual floor set by the final learning rate.

(a) 100 10−2

|z[P ]|

10−4

Sharpening sparsity with learning-rate annealing

Because the activation floor scales as O(αt ), annealing the learning rate progressively separates deactivated branches from active ones. This is shown directly in Fig. 10(a) for the XY classifier. The suppressed branches track the learning-rate scale and decrease whenever the learning rate is reduced, consistent with the scaling of Eq. (D9). When the learning-rate schedule is disabled, Fig. 10(b), the learning rate remains fixed at α = 10−2 , and the branch activation floor correspondingly stays at the same O(10−2 ) scale, overlapping with the weakest active branches. The sparse bottleneck then becomes much harder to identify. Learning-rate annealing is therefore essential in practice for obtaining a clear separation between active and deactivated branches. We compared cosine annealing, exponential decay, onecycle scheduling, and ReduceLROnPlateau using the same training budget. All four produced a clear separation between active and deactivated branches by the end of training, suggesting that the key ingredient is annealing itself rather than the schedule’s precise shape. Throughout this work, we use ReduceLROnPlateau, halving the learning rate whenever the validation metric does not improve for five epochs.

Appendix E: Unsupervised detection of transitions in the experimental sweeps: method comparison

To define the supervised classification labels used by TetrisCNN, we first estimate the transition location in each experimental dataset using five established unsupervised approaches: learning by confusion (LBC) [13], the prediction-divergence method (PDM) [15, 46, 85], principal component analysis (PCA) [81], Uniform Manifold Approximation and Projection (UMAP) [122], and distance learning (DL) [20]. For consistency, all transition estimates are reported as the interval between the two acquisition times flanking the detected change. For LBC, this interval is additionally enlarged when the inferred transition location varies across random initializations or sparsity strengths. The resulting estimates are summarized in Fig. 4(a,d).

Learning by confusion

Learning by confusion (LBC) [13] locates a transition by repeatedly imposing artificial binary labels on the

10−6 10−8 0

100

150

Epoch

(b) 100 10−1 10−2

|z[P ]|

10−3

■ ■■ ■ ■ ■□ □■ □■ ■□ ■■ ■□ ■□ ■■ □■ ■■ ■■ □■ ■■ ■■ LR

10−4 10−5 0

100

Epoch

FIG. 10: Learning-rate annealing separates active and deactivated branches. Training history of the XY classifier [rotation unconstrained variant of panel (e) of Fig. 4] with all settings kept fixed except for the learning-rate schedule. (a) With ReduceLROnPlateau, the deactivated branches follow the learning-rate scale and decrease whenever it is reduced, consistent with the O(αt ) scaling of Eq. (D10). (b) With the learning rate held fixed at α = 10−2 , the residual floor remains at the same order of magnitude and overlaps with the weakest active branches. data. For a trial transition point y ′ , snapshots are labeled according to whether their parameter lies below or above y ′ , and a classifier is trained to distinguish the two groups. Sweeping y ′ across the experimental range gives the characteristic confusion curve, whose interior accuracy maximum estimates the transition, yc ≈ arg max Acc(y ′ ). ′ y

(E1)

We use TetrisCNN itself as the classifier for each trial split. Figure 11 shows the resulting training and validation accuracies. For the XY dataset, the interior maximum is well defined and consistently falls on the third acquisition-time partition, corresponding to t ≈ 0.625 µs and bracketed by [0.5, 0.75] µs. For the Ising dataset, the maximum is broad and shallow: its location shifts between neighboring partitions depending on the random initialization, sparsity strength λmax , and whether training or validation

(a) (a)

1.00

0.9

0.98

0.7

0.9

I(y)

Accuracy

1.0

1.1

0.8 1

6 0

ŷ

t (µs)

(b)

−5

ŷ = y PT = -2.17

0.9 −7.5

0.8 0

t (µs)

FIG. 11: Learning-by-confusion transition estimate. Training (red) and validation (black) classification accuracy as a function of the hypothesized split time for (a) the Ising and (b) the XY datasets. The red dashed line marks the interior accuracy maximum of Eq. (E1). Accuracies are averaged over five random initializations; most error bars are smaller than the markers. For XY, the maximum gives t ≈ 0.625 µs, bracketed by [0.5, 0.75] µs. For Ising, the maximum is broad and shifts between neighboring partitions, yielding t ≈ 0.9 µs within [0.8, 1.2] µs. Because the experimental sweeps start close to the central confusion maximum, only the second half of the characteristic “W” curve is resolved.

accuracy is used. Across eleven values of λmax and five random initializations per value, we therefore report the median maximum, t ≈ 0.9 µs, and use the envelope of the inferred transition locations to obtain the broader interval [0.8, 1.2] µs. 2.

Prediction-divergence method

The prediction-divergence method (PDM) [15, 46, 85] trains a single regressor to predict an experimental tuning parameter y from a snapshot x, ŷ = fθ (x).

(E2)

The transition is then identified from the response of the average prediction to the true tuning parameter, I(y) =

∂⟨ŷ⟩y , ∂y

I(y)

(b) 1.6 0.8 0.0

−5.0

−2.5

0.0

2.5

Ground Truth y = δ (MHz · 2π)

15 10

ŷ

Accuracy

1.5 1.0 0.5

5 ŷ = y PT = 4.16

0 0

Ground Truth y = δXY (MHz · 2π)

FIG. 12: Prediction-divergence transition estimate. Phase indicator I(y) (top) and average prediction ⟨ŷ⟩y against the true tuning parameter (bottom) for (a) the Ising dataset, regressing δ, and (b) the XY dataset, regressing δXY . The red dashed line marks the maximum of the first-order finite-difference indicator. The corresponding transition estimates are t ≈ 1.1 µs, bracketed by [1.0, 1.2] µs (Ising), and t ≈ 0.625 µs, bracketed by [0.5, 0.75] µs (XY). In experimental parameter space, the maxima occur at δ/2π ≈ −2.2 MHz and δXY /2π ≈ 4.2 MHz, respectively.

is the divergence [85]. We use TetrisCNN as the regressor and train it with mean-squared error, with regression quality evaluated using the aggregated R2 defined in App. G. As shown in Fig. 12, the indicator peaks at the third interval, corresponding to t ≈ 1.1 µs for the Ising dataset and to t ≈ 0.625 µs for the XY dataset. For Ising we use regression of δ; for XY we regress δXY . We evaluate Eq. (E3) using a first-order finite difference, so the maximum naturally lies between two neighboring acquisition times. These two times define the uncertainty interval, giving [1.0, 1.2] µs for Ising and [0.5, 0.75] µs for XY.

(E3)

where ⟨ŷ⟩y denotes the average network prediction over all snapshots acquired at the same true parameter value y. The maximum of I(y) marks the point at which the learned prediction changes most rapidly with y. For multivariate tuning parameters, the corresponding quantity

Dimensionality reduction and clustering

As a complementary approach, we project the raw experimental snapshots into a low-dimensional representation and cluster the resulting embeddings. We use

(b)

(c)

(d) 6

(a)

PC 1

t (µs)

UMAP 2

PC 2

UMAP 2

PC 3

2 1

PC 2 UMAP 1

PC 1

UMAP 1

FIG. 13: Low-dimensional embeddings of the experimental snapshots. (a,b) Ising and (c,d) XY snapshots projected with PCA (a,c) and UMAP (b,d), colored by acquisition time. In all cases, the earliest acquisition times form a distinct group that is isolated by k-means clustering. The later-time snapshots overlap substantially and their cluster assignments no longer follow the sweep monotonically.

TABLE IV: Clustering of low-dimensional snapshot embeddings. For each dataset and dimensionality-reduction method, we report the smallest embedding and number of k-means clusters for which the initial-time cluster is isolated, together with the silhouette score and the transition interval inferred from the majority-cluster assignment. UMAP is run with n_neighbors = 250, min_dist = 0, and the cosine metric. Dataset

Reducer Components, k Silhouette Transition interval

Ising UMAP Ising PCA XY (X, Z) UMAP XY (X, Z) PCA

2, k = 3 3, k = 4 2, k = 3 2, k = 3

0.93 0.48 0.62 0.51

1.0–1.2 µs 1.0–1.2 µs 0.5–0.75 µs 0.5–0.75 µs

PCA [81], which projects the data onto directions of maximal variance, and UMAP [122], which constructs a nonlinear embedding that approximately preserves local neighborhood structure. The projected snapshots are clustered with k-means, and each acquisition time is assigned the majority cluster of its snapshots. We identify the transition where this majority assignment first changes. Figure 13 shows the resulting embeddings. For both datasets, PCA and UMAP isolate the earliest acquisition times from the remainder of the sweep. In the Ising data, the initial cluster contains the majority of snapshots up to t = 1.0 µs but not at 1.2 µs, giving the interval [1.0, 1.2] µs. In the XY data, the corresponding change occurs between 0.5 µs and 0.75 µs. At later times, the cluster assignments no longer follow the sweep monotonically, so we use the embeddings only to identify this primary change. UMAP resolves the early-time cluster in two dimensions for both datasets. Two principal components are also sufficient for XY, whereas Ising requires three principal components before the early-time cluster is separated cleanly. The corresponding silhouette scores and clustering settings are listed in Tab. IV.

Distance learning

The distance-learning (DL) method of Malyshev et al. [20] compares the full snapshot distributions at different parameter values rather than embedding individual snapshots. A discriminator network is trained to estimate pairwise squared Hellinger f -divergences between the empirical snapshot distributions. The resulting divergence matrix is then treated as a distance matrix and clustered with HDBSCAN [123]; a transition is identified at the boundary between the resulting clusters. We apply the implementation accompanying Ref. [20] to the experimental datasets. In contrast to the dense parameter grids used in the original numerical benchmarks, the experimental sweeps contain only a small number of discrete acquisition times. We therefore set the HDBSCAN minimum cluster size to 2, the smallest value that still yields meaningful clusters in this setting. As shown in Fig. 14, DL isolates the earliest acquisition points from the remainder of each sweep. The resulting boundary lies at t ≈ 1.1 µs for Ising and t ≈ 0.625 µs for XY. As for the dimensionality-reduction methods, the uncertainty is determined by the discrete acquisition grid and spans the two times flanking the cluster boundary: [1.0, 1.2] µs for Ising and [0.5, 0.75] µs for XY.

Comparison of transition estimates

All five estimates identify the same primary change in each dataset. For the Ising sweep, PDM, PCA, UMAP, and DL place the transition at t ≈ 1.1 µs, bracketed by [1.0, 1.2] µs. LBC gives a compatible but broader estimate centered at t ≈ 0.9 µs with an interval [0.8, 1.2] µs, reflecting the flat confusion maximum and its variation across fits. For the XY sweep, all five methods agree on t ≈ 0.625 µs, bracketed by [0.5, 0.75] µs. This agreement defines the locations used to construct

0.60

00 6.

00 3.

(b)

00 5.

20 4.

t (µs)

40 2.

40 1.

t (µs)

(a)

0.00

1.4

0.50 1.2

1.00 1.50

1.0

2.00

t (µs)

2.40 3.20

0.8 3.00 0.6

4.20

0.4

5.00

0.2 6.00

(c)

(d)

1.0

Confidence

6.00

Confidence

Divergence

1.40

0.5 cluster 0 cluster 1

0.0 1

0.0

0.5 cluster 0 cluster 1

0.0 6

t (µs)

FIG. 14: Distance-learning transition estimate. (a,b) Pairwise squared Hellinger f -divergence matrices between the snapshot distributions at different acquisition times and (c,d) the corresponding HDBSCAN cluster assignments and confidence for the Ising (a,c) and XY (b,d) datasets. The earliest acquisition points form a distinct cluster, yielding transition estimates of t ≈ 1.1 µs within [1.0, 1.2] µs for Ising and t ≈ 0.625 µs within [0.5, 0.75] µs for XY. Figure generated using code adapted from the repository accompanying Ref. [20].

the supervised classification labels in the main TetrisCNN analysis. For the Ising dataset, the detected change corresponds to the loss of the initially strong Z-polarization rather than to the later emergence of antiferromagnetic order, as discussed in Sec. IV. For the XY dataset, all methods consistently identify the same change near t ≈ 0.625 µs.

Appendix F: Supporting material 1.

Robustness of the learned descriptors to regularization

A key component of TetrisCNN is the L1 regularization applied to the bottleneck activations, X LL1 = λk |zk |, (F1) k

which encourages the network to identify a sparse set of descriptors sufficient for solving the classification task. Specifically, λk interpolates log-linearly between

λmin and λmax according to the pattern cardinality |Pk | [Eq. (11)]. To investigate the effect of the regularization strength, we vary the maximum penalty λmax while keeping λmin = 10−3 and the interpolation fixed. This hyperparameter controls the simplicity of the learned representation: increasing λmax penalizes larger-filter branches more strongly and therefore favors solutions based on fewer and lower-order correlators. Figure 15 summarizes the resulting branch activations and classification accuracies, averaged over ten random initializations. For the Ising dataset, Fig. 15(a), the same 1 × 1 branch is consistently selected across a broad range of λmax , once the regularization exceeds a threshold around λmax ≈ 10−1 . Then the network largely suppresses all remaining branches and converges to the one-dimensional representation discussed in the main text. Importantly, this simplification has almost no effect on predictive performance: the classification accuracy remains essentially constant throughout the entire range of sparsity strengths considered, Fig. 15(c). The value λmax = 103 used in the figures in the main text lies within this stable regime. For the XY dataset, we see a similar behavior. After reaching the penalty threshold at λmax ≈ 10−1 , the

z [P ]norm

(a)

(b)

1.0

0.5

0.0

0.936

0.99150 0.99125

(d)

0.934 −3

−2

−1

−3

log10 λmax

−2

−1

log10 λmax

FIG. 15: Robustness of branch selection to the sparsity regularization strength λmax . (a,b) Normalized bottleneck activations and (c,d) validation classification accuracy as a function of the maximum sparsity coefficient λmax for the Ising (a,c) and XY (b,d) datasets. Results are averaged over ten random initializations.

z [P ]norm

(a)

(b)

1.0

0.5

0.0

(d)

0.99175 0.99150 0.99125

0.945 0.940 −3

−2

−1

log10 λmax

−3

−2

−1

log10 λmax

FIG. 16: Robustness of branch selection to the sparsity regularization strength λmax for filters that are not rotationally invariant. (a,b) Normalized bottleneck activations and (c,d) validation classification accuracy as a function of the maximum sparsity coefficient λmax for the Ising (a,c) and XY (b,d) datasets. Results are averaged over ten random initializations.

network consistently suppresses the redundant branches while retaining the same three dominant descriptors identified in the main text. Once again, this reduction in representation complexity causes only a negligible change in classification accuracy, as seen in Fig. 15(d), and the main-text value of λmax = 103 lies within this stable region. Taken together, these results demonstrate that the physically interpretable descriptors identified by TetrisCNN are not sensitive to the precise choice of sparsity strength. Instead, the sparsity regularization acts primarily as a model-selection mechanism, selecting a simple representative from a large family of functionally equivalent solutions with nearly identical predictive performance. This observation is consistent with the Rashomon effect [42]: many distinct internal representations achieve comparable classification accuracy, while the sparsity penalty biases the network toward the simplest among them.

Searching for the relevant scale

The main text restricts the branch dictionary to patterns no larger than 2 × 2 (Tab. I). Before committing to this choice, we tested whether TetrisCNN would favor a larger spatial scale if given the option. We repeat the same classification setup and λmax sweep used to fix the label boundary (App. E 1), but replace the branch dictionary with four deliberately coarse-grained kernels: 1 × 1, 2 × 2, 4 × 4, and the full system size (8 × 8 for the Ising data, 6 × 7 for the XY data). Figure 17 shows the resulting branch activations and validation accuracy as a function of λmax , averaged over ten random initializations. For both datasets, the network consistently suppresses kernels larger than 2 × 2 once the sparsity penalty exceeds a modest threshold with little effect on classification accuracy. For the XY data, suppressing the 4 × 4 branch is accompanied by a small decrease in classification accuracy of approximately 1%, suggesting this branch captures some

z [P ]norm

(a)

(b)

1.0

0.5

0.0

(c)

(d)

0.991 0.990

Acc

0.94 0.93 −3

−2

−1

log10 λmax

−3

−2

−1

log10 λmax

FIG. 17: Searching for the relevant scale for the Ising (first column) and XY model (second column). (a,b) Normalized bottleneck activations and (c,d) validation classification accuracy as a function of the maximum sparsity coefficient λmax for the Ising (a,c) and XY (b,d) datasets. Averaged over 10 random initializations.

5 2.

75 1.

25 1.

Interestingly, the best unconstrained network trained on the XY dataset in Fig. 16(b),(d) achieves approximately one percentage point higher validation accuracy than the best rotationally invariant network. As shown in Fig. 18, this improvement is concentrated near the phase transition, indicating that the experimental snapshots contain some additional classification-relevant information that is not rotationally invariant. Nevertheless, the small accuracy difference and agreement in the selected pattern types support our use of the rotationally invariant model for the physical interpretation.

unconstrained invariant (C4 )

To test the effect of imposing the C4 rotational symmetry of the square lattice, we compare the rotationally invariant TetrisCNN used in the main text with an unconstrained variant across the same sweep of λmax . As shown in Fig. 16, the unconstrained variant identifies the same types of relevant branches as the rotationally invariant network (Fig. 15). For Ising, the learned representations are identical. For XY, the rotationally invariant model retains three branches and the unconstrained model retains five, corresponding to separate rotational variants of the same patterns. We therefore focus on the rotationally invariant network in the main text, as it provides an equivalent but lower-dimensional and simpler representation.

The X- and Z-basis measurements in the XY dataset correspond to independent experimental realizations. Consequently, pairing individual X- and Z-basis snapshots into two-channel inputs is arbitrary. To verify that the interpretation reported in the main text does not depend on this choice, we repeat the complete TetrisCNN analysis for ten independent random re-pairings of the two measurement bases, using the same hyperparameters and sparsity strength as in the main-text analysis. For each re-pairing, we train TetrisCNN from scratch, identify the surviving branches and their leading Boolean-Fourier components, and independently fit the quadratic decision boundary in the resulting bottleneck space, as described in App. H 1. Because multiplying a decision-boundary equation by a nonzero constant leaves the boundary un-

Rotationally invariant and unconstrained TetrisCNN

Robustness to random pairings of X and Z snapshots

Error rate (%)

useful nonlocal information. Nevertheless, the vast majority of the predictive signal comes from the local branches. These results indicate that local patterns are sufficient to capture nearly all information relevant to the classification tasks considered here, justifying our restriction to patterns no larger than 2 × 2 in the main text.

t (µs)

FIG. 18: Classification validation error across the XY sweep for rotationally invariant and unconstrained TetrisCNN. The unconstrained network achieves approximately 1% lower validation error, with the improvement concentrated near the phase transition around 0.625 µs.

(a)

Ising

5. 8

5 5. 4

4. 6

4. 2

3. 8

3 3. 4

2. 6

2. 2

1. 8

0. 6

1 1. 4

Error rate (%)

t (µs) (b)

XY, XZ

20 10

5 1. 5 1. 75

1. 2

5 0. 7

0. 2

t (µs) Error rate (%)

changed, we normalize the fitted equations to the coeffiZ cient of Crot [■■] before comparing their coefficients. We then report the mean and standard deviation of each normalized coefficient across the ten re-pairings in Tab. V, compare them with the main-text fit in Eq. (14), and quantify the deviation of each main-text coefficient from the re-pairing mean in units of the corresponding standard deviation. The resulting bottleneck representations are highly stable, as presented in Tab. V. Across all ten re-pairings, TetrisCNN selects the same three branches, whose leading Boolean-Fourier components correspond to C X [■], Z Z ■□ Crot [■■], and Crot [□■]. Moreover, the fitted quadratic decision boundary has the same functional form across runs. As shown in Tab. V, after normalization, all coefficients of the main-text decision boundary lie within 1.2 standard deviations of their means across the ten re-pairings, which demonstrates that both the selected physical correlators and the decision boundaries reported in the main text are robust to the arbitrary pairing of Xand Z-basis snapshots.

Error rate (%)

(c)

XY, X

XY, Z

XY, XZ

To quantify the information available in the two measurement bases, we construct XY datasets containing the same total number of experimental snapshots, n = 14471, sampled at the same sweep points: n X-basis snapshots, n Z-basis snapshots, and n/2 snapshots from each basis in the combined XZ dataset. We train TetrisCNN independently on each dataset using five random initializations. As shown in Fig. 19, the Z-basis snapshots alone contain most of the information required for classification, reaching 91.19% validation accuracy, compared with 70.28% for the X basis alone. Combining the two bases further increases the accuracy to 93.07% (with standard deviations across five random initializations below 0.01% for all three networks), despite using only half as many paired inputs, with the largest improvement occurring near the identified phase transition. This suggests that the X-basis snapshots contain complementary informaZ Z ■□ C X [■]2 C X [■] Crot [■■] Crot [□■] const 0.590 0.278 1 −0.761 0.347 Coefficients ±0.092 ±0.060 ±0.060 ±0.018 Main-text 0.680 0.206 1 −0.815 0.346 Deviation [σ] 0.98 −1.2 −0.9 −0.06

Terms

TABLE V: Robustness of the quadratic XY decision boundary to X-Z snapshot pairing. Mean and standard deviation of the normalized coefficients obtained from ten independently re-paired and retrained networks, compared with the main-text result in Eq. (14). All main-text coefficients lie within 1.2σ of the re-pairing ensemble.

5 2.

1 1. 25 1. 5 1. 75

0. 25 0. 5 0. 75

Contribution of the X and Z measurement bases to classification

t (µs)

FIG. 19: Classification validation error rates per time for TetrisCNN (a) from the main text, trained on the Ising dataset, (b) from the main text, trained on the XY dataset, and (c) trained on the XY dataset when using only X-snapshots, only Z-snapshots, and both XZ snapshots in the comparable setup with the same total number of snapshots. All networks make their biggest mistakes around the detected crossover (Ising) or phase transition (XY), partly because single-snapshot classification limits accuracy. In the XY dataset in panel (c), the network trained on Z snapshots only has better accuracy than the one trained on X snapshots only (91.19% and 70.28%, respectively). Information from the X basis helps mostly around the phase transition, allowing the network trained on both bases to increase its total validation accuracy to 93.07%.

tion about the transition, consistent with their sensitivity to the XY order parameter discussed in the main text. To understand what information is extracted from each basis, we additionally fit decision boundaries to the TetrisCNNs trained on the X- and Z-basis snapshots separately. Across all five initializations, the X-only classifier learns a quadratic boundary in the X magnetization, (C X [■])2 + 0.85 C X [■] − 0.36 = 0,

(F2)

consistent with the relation between (C X [■])2 and the XY order parameter discussed in the main text. However, its substantially lower classification accuracy shows that this information alone is insufficient to reliably distinguish

29 ■□ □■

Classification

1.0 0.5 0.0 2

■■ ■□

■□ ■■

□■ ■■

■■ □■

Regression: δ

(g)

0 −5

t (µs)

(b)

□■ ■□

(d) δ/2π (MHz)

softmax[logit0 ]

(a)

■ ■

■■

Ω/2π (MHz)

■

t (µs)

(h)

0.1

0.5

z[P ]

Regression: Ω

t (µs)

(e)

■■ ■■

0.0

0.0 2

(c)

t (µs)

(f)

(i)

10−2 10−6

100

|z[P ]|

10−6

t (µs)

10−3

|z[P ]|

10−2

t (µs)

10−6

Epoch

100

Epoch

150

100

Epoch

FIG. 20: Task dependence of the learned TetrisCNN representation. Comparison of classification (a–c), regression of δ (d–f), and regression of Ω (g–i) for the Ising dataset, using the same TetrisCNN architecture and branch regularization. Top row: classification probability or regression target versus sweep time. Middle row: branch activations z[P ] versus sweep time. Bottom row: mean absolute branch activations during training.

the two classes from individual snapshots. In contrast, the Z-only classifier learns the stable linear boundary Z Z ■□ 0.776 Crot [■■] − Crot [□■] + 0.046 = 0,

(F3)

involving the same two Z-spin correlators that dominate the combined XZ classifier. Thus, TetrisCNN supplements the physically expected information from the squared X-magnetization with Z-basis correlators that provide substantially stronger discrimination at the level of individual snapshots. 6.

Task dependence of the learned descriptors

The main-text Discussion argues that the correlators selected by TetrisCNN depend on the learning objective. Here we demonstrate this directly, for the Ising data and rotationally unconstrained network, by comparing the branches selected for classification with those selected by the regression (implemented within the predictiondivergence method, cf. App. E 2) of the two experimental tuning parameters separately. Figure 20 shows classification (first column, as in the main text), regression of

δ (second column), and regression of Ω (third column), using otherwise equivalent TetrisCNN architectures and regularization (λmax = 103 ). Classification and regression of δ both predominantly select the single-site branch z[■]. In contrast, when regressing Ω, z[■] is suppressed and the network instead □■ relies on z[■■], z[■ ■] and z[■□] branches. The representation selected by the network therefore reflects not only the underlying data but also which information is required by the learning objective. In particular, classification with TetrisCNN favors the simplest correlators sufficient to distinguish the two classes, whereas regression can favor different correlators informative of continuous variation of the target along the experimental sweep.

Complete mapping from filters to correlators

Figure 21 presented in this appendix reproduces Fig. 5 of the main text without pruning: every subpattern correlator C[P ′ ], P ′ ⊆ P enters the multilinear fit, whereas the main-text version collapses the smaller positive and negative contributions into the two “other” categories for

30 z [] ≈ −1.56 C

(a)

X: Z:

+ 0.02 | R2 = 1.00

−1 ∅

X: Z:

z [] ≈ −0.21 Crot

(b)

X: Z:

− 0.03 | R2 = 0.97

0.0 −0.1 −0.2

∅

X: X: X: X: X: X: X: X: X: X: X: X: X: X: X: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z:

(c)

≈ 0.16 Crot

X: Z:

− 0.01 | R2 = 0.98

0.1

0.0 ∅

X: X: X: X: X: X: X: X: X: X: X: X: X: X: X: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z: Z:

FIG. 21: Mapping from filters to correlators via multilinear regression for three active branches in the XY model, enabled by the Boolean Fourier expansion, Eq. (7). We only consider correlators C[P ′ ], P ′ ⊆ P defined by patterns P ′ □ that are contained in the filter pattern P (see Tab. I). Note that we need to take into account both C[■ □] and C[■] because edge effects make these two correlators have different values. This is the full version of the Fig. 5 in the main text.

readability. The two fits agree by construction, since the main-text figure is obtained from this full decomposition by pruning; it is reproduced here in full for transparency.

Appendix G: R-squared for multivariate and aggregate output

When assessing the quality of a regression model where the predicted output is multivariate (e.g. when predicting both δ and Ω for the TFI model in the main text), common libraries (PyTorch [124] and scikit-learn [125]) allow a choice between “uniform weighting” and “variance weighting” the R2 of the two outputs. We argue here that variance weighting is appropriate if the canonical interpretation of R2 is to be preserved (i.e. as a proportion of variance in the data explained by the model).

That is, for multivariate y ∈ RK , with N datapoints for each k ∈ [K], R2 for the model as a whole is defined as: 2 P P (n) (n) ŷk − yk k n∈Γ SS k res 2 =1− P P Rtot =1− 2 SStot (n) ŷ − ȳ k k n∈Γk k (G1) P where ȳk = N1k n∈Γk yn is the mean for group k. We use the subscript “tot” to distinguish from R2 for an individual output yk : 2 P (n) (n) ŷ − y k n∈Γ k k SSres k Rk2 = 1 − =1− P 2 k SStot (n) n∈Γk yk − y k k Using SSres = SStot 1 − Rk2 one finds after some alge-

31 bra that: mean |coef|

40 20 0

1.00

where P the sums run over2 unique values of y and ỹ = 1 y y is their mean. Ragg quantifies the variance in y N explained by groups of snapshots sharing a label, whereas the standard R2 measures the variance explained by individual datapoints. For multivariate y it is variance weighted, following the argument above. Appendix H: Fitting the decision boundary to the classifier predictions in the bottleneck space

C (inverse L1 penalty)

FIG. 22: Regularization path for the polynomial surrogate of the decision boundary of a single TetrisCNN. Coefficients of the polynomial features are shown as a function of the inverse L1 regularization strength C, together with the five-fold cross-validation accuracy. Four terms dominate as the regularization is relaxed: z[■]2 , z[■], z[■■], and z[■□ □■]. Error bars show the standard deviation over five cross-validation realizations. where σ(x) = 1+e1−x is the sigmoid function. Since σ(0) = 1/2, the analytical approximation to the decision boundary is the zero-level set F (z) = β0 + β T ϕ(z) = 0.

Fitting the decision boundary

For the XY dataset, where the TetrisCNN decision boundary is a surface in the 3D bottleneck space rather than a single threshold as for Ising, extracting an analytical approximation requires choosing its functional form and the terms included in it. To make these choices systematic, we construct a logistic-regression surrogate of the trained TetrisCNN classifier in its interpretable bottleneck space. This approach is conceptually related to surrogatemodel explainability methods such as LIME [126], but here we fit the surrogate globally rather than locally and use it to extract an analytical approximation of the decision boundary in an already interpretable latent space. For each snapshot S, the surrogate takes the three active bottleneck coordinates z(S) = z[■], z[■■], z[■□ (H1) □■] and their polynomial combinations up to degree two, while its target is the class ŷTetris (S) predicted by TetrisCNN rather than the ground-truth label. Defining ϕ(z) = (z1 , z2 , z3 , z12 , z22 , z32 , z1 z2 , z1 z3 , z2 z3 ),

(H2)

we fit the logistic-regression surrogate h i p(ŷTetris = 1 | z) = σ β0 + β T ϕ(z) ≡ σ[F (z)],

(H3)

st 50 ra 0 in ed on

un c

The same distinction between individual and grouped datapoints arises for the aggregated R2 tracked alongside the MSE loss in the prediction-divergence method (App. E 2), P 2 y (y − ⟨ŷ⟩y ) 2 (G2) Ragg = 1 − P 2 , y (y − ỹ)

5 100 0

±1 std over 5 regression seeds

0.95

0. 00

z22 z2 z3 z32

01 0. 00 0. 05 00 1 0. 00 0. 5 01 0. 0 0.5 1

1 X 2 2 Rk ̸= Rtot . K

z12 z1 z2 z1 z3

5 10

I.e. one obtains a ‘variance-weighted’ average of the in2 2 dividual = P 2 Rk . I.e. it is clear that the “average” R 1 2 R is in general different from the “total” R , dek k K fined by (G1):

z1 z2 z3

0. 5 1

X SS k SS k tot P res k′ = Rk2 . SS SS ′ tot tot k k

CV acc.

2 Rtot =1−

(H4)

Accordingly, throughout the main text and below, we report the polynomial decision function F (z) rather than the logistic probability itself. All model selection is performed using only the training snapshots, with five-fold stratified cross-validation. The reported accuracy is evaluated on the held-out validation snapshots and measures agreement with the TetrisCNN predictions. To identify the polynomial terms most relevant to the decision boundary, we examine an L1 -regularized logistic-regression path, varying the inverse regularization strength C from strong regularization toward the effectively unregularized limit. As shown in Fig. 22, four terms dominate: z[■]2 , z[■], z[■■], and z[■□ □■]. Through the Boolean-Fourier mapping (Fig. 5), these corZ respond predominantly to (C X [■])2 , C X [■], Crot [■■], Z ■□ and Crot [□■], respectively. The same four terms are consistently recovered across independent TetrisCNN initializations (Fig. 23), motivating the linear and quadratic approximations considered below. Further fitting details are given in Tab. VI. 2.

Decision boundary: plane vs quadratic surface

Having established the dominant terms, we next ask how much of the TetrisCNN decision rule progressively

st 50 ra 0 in ed

0.00

t (µs)

un c

5 100 0

0. 5 1

0. 00

5 10

±1 std over 5 model inits

0.95

0.01

1.00

01 0. 00 0. 05 00 1 0. 00 0. 5 01 0. 0 0.5 1

CV acc.

0.02

Plane Square

0.03

0 1. 25

z22 z2 z3 z32

z12 z1 z2 z1 z3

0. 5 0. 75

z1 z2 z3

0. 0 0. 25

mean |coef|

Error rate (1 - accuracy)

C (inverse L1 penalty)

FIG. 23: Stability of the regularization path across several TetrisCNN run initializations. Regularization paths averaged over five independently initialized TetrisCNN models. The same four dominant polynomial terms are consistently recovered, demonstrating that the inferred functional form of the decision boundary is stable across TetrisCNN random initializations. TABLE VI: Logistic-regression settings used to fit the polynomial decision-boundary approximations. Setting

Value

Target TetrisCNN predicted class Preprocessing StandardScaler Cross-validation 5-fold stratified Final regularization L2 , C = 104 Solver lbfgs Class weight balanced Max. iterations 2 × 104 Evaluation held-out validation set

near the identified phase transition, where the quadratic approximation reduces the disagreement with the classifier from approximately 3% for the linear boundary to zero, as shown in Fig. 24. This shows that the two Z-correlators provide the dominant classification signal, while the nonlinear contribution from C X [■] refines the decision boundary, particularly in the vicinity of the phase transition.

Appendix I: Limitations of symbolic regression

more expressive physical approximations can capture. Remarkably, the linear approximation using two Zcorrelators alone already reproduces the classifier predictions with 96.53% accuracy: Z Z ■□ Crot [■■] − 0.879Crot [□■] + 0.41 = 0 .

FIG. 24: Linear and quadratic approximations to the XY decision boundary. Disagreement with the TetrisCNN predictions as a function of sweep time for the linear and quadratic approximations to the decision boundary in Fig. 6(b). The improvement obtained by including (C X [■])2 is concentrated near the identified transition at approximately 0.625 µs.

In the original TetrisCNN pipeline [69], we used symbolic regression to obtain an analytical approximation of the task network acting on the interpretable bottleneck. Here we revisit this approach and find that fitting the full task-network output can yield symbolic expressions that agree well with network predictions in distribution but extrapolate poorly.

Search space and evaluation

(H5)

Thus, the majority of the classifier’s decision rule can be explained without the X-magnetization contribution. Including C X [■], the full linear approximation becomes

We use PySR [127], which searches over symbolic expressions using an evolutionary algorithm. We restrict the operator basis to

Z Z ■□ 0.22C X [■] + Crot [■■] − 0.8Crot [□■] + 0.39 = 0 , (H6)

{+, ×, negation, (·)2 , (·)3 , (·)4 },

and reproduces the classifier predictions with 96.85% accuracy. Finally, only allowing a quadratic dependence on C X [■] as in the main text in Eq. (14), improves the agreement more significantly. While the overall improvement seems modest (accuracy increases to 98.30%), it is concentrated

favoring low-order polynomial forms natural for local spin observables. Multiplication is restricted to multiplication by a constant, powers cannot be nested, and expression complexity is explicitly penalized. We use mean absolute error as the fitting loss and a parsimony penalty of λ = 0.01. Expanding the operator basis, for example by

(I1)

33 including division or exponentials, does not qualitatively change the conclusions below. Each symbolic expression is fitted to reproduce the output of the trained task network from the active bottleneck variables: the classification logit for the Ising classifier or the predicted transverse field Ω for the Ising regression 2 task. In-distribution agreement is quantified by Rval on validation snapshots. To test whether the symbolic form captures the task network beyond the training distribution, we evaluate it on synthetic spin configurations not encountered during training. Because the targets for these configurations are provided by the output of the trained network itself, we are free to probe arbitrary spin configurations. We consider paramagnetic, ferromagnetic, and antiferromagnetic families, each generated from a fixed representative configuration with four random spin flips. These configurations are passed through the original network and the symbolic surrogate, and the agreement is evaluated separately for each family. Alongside R2 , we report the Spearman correlation ρ, which distinguishes quantitatively poor fits that nevertheless preserve the ordering of the network predictions.

Classification: full network output versus decision boundary

The Ising classifier provides the simplest setting, since its bottleneck collapses to the single branch z[■]. Symbolic regression is trained to reproduce the classifier logit over the full bottleneck domain. Several expressions on the PySR Pareto front achieve essentially identical classification accuracy, while their agreement with the logit ranges 2 from Rval = 0.50 to 1.00 and their out-of-distribution behavior differs substantially, with the strongest degradation 2 occurring for the ferromagnetic probe family (Rval = 0.56– 0.85). Importantly, this is a stricter requirement than reproducing the classifier itself, as we do in the main text and describe in detail in App. H 1. A classification surrogate only needs to recover the location of the decision boundary, and the detailed variation of the logit away from that boundary is irrelevant to the predicted class. We therefore evaluate the polynomial decision-boundary approximation used in the main text on the same out-of-distribution configurations. It reproduces the TetrisCNN class predictions with 100% accuracy on all three out-of-distribution probe families. The decision-boundary fit generalizes perfectly to the out-of-distribution probes, whereas the symbolic-

regression fit to the full network logit does not. This may be because reproducing the full logit promotes more complicated, less generalizable functions, or because symbolic regression additionally depends on the user-defined operator grammar and complexity prior. What is clear, however, is that a controlled approximation of the decision boundary provides the more robust surrogate in this setting.

Higher-dimensional bottleneck: failure to extrapolate

The limitations become clearer for the Ising Ωregression task. Here the single-site branch is suppressed (App. F 6), and the learned representation is distributed across several two-body branches. No expression on the PySR Pareto front is simultaneously simple, accurate in distribution, and robust across the out-of-distribution probes. 2 In particular, Rpara is negative for every expression in Tab. VII, and the simpler candidates also fail on the ferromagnetic and antiferromagnetic families. The Spearman correlation shows that these failures can differ qualitatively: for example, one expression reverses the ferromagnetic ordering (ρferro = −0.90), whereas another preserves it (ρferro = 0.90) despite still having poor R2 . The out-of-distribution predictions are not arbitrary [Fig. 25(d)]: they remain structured by configuration family and largely within the numerical range encountered in distribution. Nevertheless, their quantitative disagreement with the network demonstrates that a good symbolic fit on the available validation data does not identify a unique or robust functional form. Taken together, these results motivate a more restricted role for symbolic regression in TetrisCNN. For classification, fitting the full task-network output is unnecessarily demanding when the scientific object of interest is the decision boundary: a simple polynomial fit to that boundary generalizes perfectly across the out-of-distribution probes considered here. For regression, where no analogous boundary reduction is available, symbolic regression can fit the network well in distribution yet extrapolate poorly. Moreover, the selected symbolic form depends on the user-defined operator grammar and complexity prior. We therefore use architecture-native interpretation to identify the physical bottleneck variables, and, when possible, fit only the simplest task-relevant object in that interpretable space rather than attempting to reconstruct the entire task network symbolically.

(b) Classifier, out-of-distribution

0.4

−20

0.3

−50

−40

−30

−20

−10

Network prediction fNN

0.2

Paramagnetic Antiferromagnetic Ferromagnetic

−40

−50

−40

−30

−20

−10

0.1 10

Network prediction fNN

0.0

(d) Ω-regression, out-of-distribution

1.0

−0.1

0.8

−0.2

Residual fSR − fNN

SR prediction fSR

(a) Classifier, in-distribution

0.6 −0.3

0.4 Paramagnetic Antiferromagnetic Ferromagnetic

0.2 0.2

0.4

0.6

0.8

1.0

Network prediction fNN

0.2

0.4

0.6

0.8

−0.4

1.0

Network prediction fNN

FIG. 25: In-distribution agreement does not guarantee out-of-distribution robustness. Symbolic-regression prediction fSR against network prediction fNN for the equations marked † in Tab. VII. Top: Ising classifier; bottom: Ising Ω regression. Left panels show validation data and right panels the paramagnetic, ferromagnetic, and antiferromagnetic probe families. Points are colored by the signed residual fSR − fNN .

TABLE VII: Symbolic-regression fits and out-of-distribution generalization. For each equation, we report the 2 validation fit to the network output, Rval , and the coefficient of determination R2 and Spearman correlation ρ on paramagnetic, ferromagnetic, and antiferromagnetic probe families. Rows within each task correspond to expressions of increasing complexity along a single PySR Pareto front. † marks the expressions shown in Fig. 25. paramagnetic ferromagnetic antiferromagnetic Task

Equation

2 Rval R2

Ising clf. Ising clf. Ising clf. Ising clf.

8.60 z[■] −z[■]2 + 6.35 z[■] − 1.22 † −1.55 z[■]2 + 5.69 z[■] − 1.19 −0.419 z[■]3 − 1.73 z[■]2 + 6.73 z[■] − 0.481

0.50 1.00 0.95 0.99 0.99 0.96 1.00 1.00

1.00 1.00 1.00 1.00

0.84 0.85 0.56 0.76

0.95 0.95 0.95 0.95

1.00 1.00 0.97 1.00

0.96 0.96 0.96 0.96

Ising reg. (Ω) 11.5 z[■ 0.82 −5.9 ■] + 0.386 3 ■] + 0.407 Ising reg. (Ω) −1.11 × 103 z[■ ] + 13.8 z[ 0.91 −7.1 ■ ■ 3 ■] + z[□■] + 0.406 † 0.92 −6.3 Ising reg. (Ω) z[■■] − 1.22 × 103 z[■ ] + 12.9 z[ ■ ■ ■□

−0.31 −0.31 −0.25

−4.8 −0.04 0.38

−0.90 0.90 0.87

−0.62 −0.26 0.12

0.72 0.72 0.76

Record · ID 978421 · SHA-256 d0c9c225c2d011c2
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.