Parallel Scan Recurrent Neural Quantum States for Scalable Variational Monte Carlo Ejaaz Merali,1, 2, ∗ Mohamed Hibat-Allah,3, 4, ∗ Mohammad Kohandel,3 Richard T. Scalettar,1 and Ehsan Khatami2, †
arXiv:2605.13807v1 [cond-mat.str-el] 13 May 2026
2
1 Department of Physics and Astronomy, University of California, Davis, California 95616, USA Department of Physics and Astronomy, San José State University, San José, California 95192, USA 3 Department of Applied Mathematics, University of Waterloo, Waterloo, ON N2L 3G1, Canada 4 Vector Institute, Toronto, Ontario, M5G 0C6, Canada (Dated: Thursday 14th May, 2026)
Neural-network quantum states have emerged as a powerful variational framework for quantum many-body systems, with recent progress often driven by massively parallel architectures such as transformers. Recurrent neural network quantum states, however, are frequently regarded as intrinsically sequential and therefore less scalable. Here we revisit this view by showing that modern recurrent architectures can support fast, accurate, and computationally accessible neural quantum state simulations. Using autoregressive recurrent wave functions together with recent advances in parallelizable recurrence, we develop variational ansätze, called parallel scan recurrent neural quantum states (PSR-NQS), which can be trained efficiently within variational Monte Carlo in one and two spatial dimensions. We demonstrate accurate benchmark results and show that, with iterative retraining, our approach reaches two-dimensional spin lattices as large as 52 × 52 while remaining in agreement with available quantum Monte Carlo data. Our results establish recurrent architectures as a practical and promising route toward scalable neural quantum state simulations with modest computational resources.
I.
INTRODUCTION
Neural-network quantum states (NQS) have emerged as a powerful framework for the variational study of quantum many-body systems, combining the expressive power of modern machine learning architectures with the flexibility of variational Monte Carlo (VMC) methods [1– 5]. Over the past few years, a broad range of neural ansätze have been explored, including convolutional networks, autoregressive models, recurrent neural networks (RNNs), and more recently transformer-based architectures [6–18]. These developments have substantially expanded the scope of NQS approaches, allowing increasingly accurate simulations of strongly correlated systems in one and two spatial dimensions [14, 16, 19–23]. Among these architectures, transformers have attracted particular attention because of their strong empirical performance and their compatibility with massively parallel modern hardware [12, 16, 17, 24]. Recent developments, however, suggest that the conventional contrast between “slow RNNs” and “fast transformers” is no longer as sharp as previously assumed. While RNNbased approaches have traditionally been viewed as intrinsically sequential and therefore comparatively slow, several recent families of recurrent or recurrence-like sequence models have shown that this limitation can be substantially mitigated. In particular, state space models (SSMs), linear recurrent units (LRUs), and simplified gated recurrent architectures [25–32] exploit structured recurrences that can be evaluated efficiently using parallel scan algorithms, thereby recovering substantial
∗ Co-first author. † [email protected]
parallelism while retaining favorable modeling properties of RNNs. Moreover, certain autoregressive transformer formulations admit an iterative state-update representation closely related to recurrence, further blurring the boundary between recurrent and attention-based sequence models [33, 34]. The recurrent feature of RNNs is especially relevant in the context of quantum many-body modeling, where architectural inductive biases can play a central role. In particular, recurrent models process configurations through repeated state updates along the chosen ordering of the lattice, which can provide an intrinsic notion of relative separation along the sequence without requiring distance information to be imposed entirely through external positional encodings [9, 35]. Although such distancedependent behavior is not guaranteed for arbitrary RNN parameterizations, common stable recurrent updates, including diagonal recurrences, can give rise to naturally decaying correlations with relative distances in recurrent neural quantum states while still allowing both shortand long-range correlations to be represented [35]. By contrast, in transformer-based approaches, distance dependence is typically introduced through additional design choices, such as positional encodings, or explicit distance-decay mechanisms [20]. This makes recurrent architectures particularly appealing for many-body problems, where distance and correlation structure are tied in the underlying quantum many-body physics. Motivated by these developments, in this work, we revisit RNN wave functions from the perspective of scalability, efficiency, and practical computational accessibility. Our goal is to show that modern recurrent architectures can support fast and accurate wave function representations in both one and two spatial dimensions, while remaining computationally lightweight. By harnessing recent advances in parallelizable recurrent models, in-
2 cluding SSM-inspired linear recurrences and simplified gated recurrent updates, we develop a parallel scan recurrent neural quantum state, dubbed PSR-NQS, that can be trained efficiently within the VMC framework. In particular, we demonstrate that these constructions can reach systems as large as 52×52 spins using less expensive computational resources than state-of-the-art methods, thereby providing a competitive alternative to computationally heavier large-scale architectures. Our aim here is not to argue that recurrent architectures are comparable in performance and accuracy with transformer-based models. Rather, we focus on the more specific question of whether modern parallelizable recurrent architectures can reduce the scalability gap, while retaining useful inductive biases that are well-suited to lattice quantum systems. More broadly, our results support the view that scalable NQS simulations need not rely exclusively on computationally intensive architectures or large hardware budgets. Instead, recurrent models offer a compelling combination of favorable inductive bias, algorithmic efficiency, and practical computational accessibility. We therefore believe that they constitute a promising route toward large-scale many-body simulations with NQS. This paper is organized as follows. In Sec. II, we first review the recurrent sequence-modeling ingredients that motivate our PSR-NQS, including SSMs, LRUs, and simplified gated recurrences. We then specify the concrete PSR-NQS implemented in this work, namely a 1D LRU wave function and a 2D minGRU wave function. In Sec. III, we present benchmark studies of the PSRNQS in terms of numerical speedup, ground-state estimation, and finite-size scaling toward the thermodynamic limit for the one-dimensional transverse-field Ising model (TFIM) and the square-lattice Heisenberg model. II. A.
METHODS
State Space Models
Recent developments in sequence modeling have shown SSMs to be effective at capturing long-range dependencies while allowing autoregressive generation to be performed sequentially with one state update per input token and without storing the full sequence history [27, 28]. SSMs are defined starting from a first-order ordinary differential equation. Given an input x(t) ∈ Rd with embedding dimension d, we map it to a hidden state h(t) ∈ Fdh with memory hidden state size dh (where F may be R or C), before finally producing an output state y(t) ∈ Rd : d h(t) = Āh(t) + B̄x(t) dt y(t) = Re{Ch(t)} + Dx(t)
(1) (2)
where Ā ∈ Fdh ×dh , B̄ ∈ Fdh ×d , C ∈ Fd×dh , and D ∈ Rd×d . This continuous-time system is then discretized,
commonly using the bilinear method [36] or the zeroorder hold method [30]. Applying the zero-order hold discretization, with time step ∆, gives the following: ht = Aht−1 + Bxt−1 (3) yt = Re{Cht } + Dxt−1 (4) with A = exp ∆Ā , B = (∆Ā)−1 (exp ∆Ā − 1)∆B̄, where the exponentials are matrix exponentials. While the matrix exponential is computationally expensive in general, it has been found that in practice Ā may be taken to be diagonal. Additionally, the matrix B may sometimes be parameterized directly [37]. Different SSM architectures vary in their parameterizations of A, B, C, with more recent architectures allowing these to depend on the input xt−1 [30]. Due to the diagonal structure of the matrix A and the linearity of the recurrence relation, Eq. (3) can be parallelized using the parallel scan algorithm [38], resulting in an overall recurrence depth of O(log N ) for a sequence of length N . This scan-based evaluation is one of the main algorithmic ingredients used in this work: whenever the recurrent update can be written as an associative affine recurrence, we replace the usual sequential loop by a parallel scan. We must emphasize that due to the autoregressive decomposition, the sampling process will still take O(N ) time. However, in VMC we often need to perform O(N ) additional forward passes of the wave function ansatz in order to compute the amplitudes of connected offdiagonal basis states of the local energy [1]. Although the additional forward passes over connected configurations are, in principle, embarrassingly parallel, our implementation does not rely on this form of parallelism. Parallelizing over all O(N ) connected configurations would require storing them simultaneously, resulting in O(N 2 ) peak memory usage per local-energy evaluation. This is prohibitive for large systems. Instead, we exploit the parallel scan structure, which parallelizes the recurrent computation itself, while sequentially iterating over the connected off-diagonal basis states. Since the local-energy evaluation is typically the dominant cost in VMC training [39], this reduces the effective scaling from the usual O(N 2 ) time to O(N log N ) time, while retaining O(N ) peak memory usage. B.
Linear Recurrent Unit
Building on this SSM perspective, one of the models we use in this study is the LRU [29], a simple architecture that nevertheless demonstrates strong long-range reasoning capabilities [29]. It specializes the discrete SSM recurrence in Eq. (3) to the case of a diagonal complex-valued propagator, whose entries are explicitly parameterized in terms of decay and phase. The LRU recurrence relation
3 is defined as ht = exp −eν + ie
θ
⊙ ht−1 + exp(γ) ⊙ (Bxt−1 ),
yt = Re{Cht } + D ⊙ xt−1 ,
(5) (6)
where ν, θ, γ ∈ Rdh are exponentiated to promote optimization stability, B ∈ Cdh ×d , C ∈ Cd×dh , D ∈ Rd , and ⊙ represents elementwise multiplication. Because the recurrence is linear and diagonal, it can be parallelized efficiently using parallel scan techniques [38]. In our implementation, this observation is used directly, where the hidden states of the 1D LRU are computed by a parallel scan rather than by iterating sequentially over the lattice sites. Note that the complex phases and exponential decay provide a natural mechanism for representing oscillatory and long-range behavior. In this sense, the LRU keeps the recurrent dynamics themselves as simple as possible, delegating any additional nonlinearity to the blocks surrounding the recurrence (see Sec. II D). C.
minGRU
Having introduced SSMs and LRUs as examples of parallelizable linear recurrences, we next consider a simplified gated recurrent architecture, the minGRU [31], which is a simplified variant of the gated recurrent unit (GRU) [40]. The minGRU recurrence is defined as: zt = σ(Wz xt−1 + bz )
(7)
h̃t = Wh xt−1 + bh
(8)
ht = zt ⊙ ht−1 + (1 − zt ) ⊙ h̃t
(9)
dh ×d
dh
where Wz , Wh ∈ R and bz , bh ∈ R . Unlike the LRU, minGRU does not use complex hidden states and does not include a built-in output map back to dimension d. Its key distinction is that the gate zt makes the recurrence input-dependent, so the hidden-state update is no longer strictly linear with fixed coefficients. This preserves a simple recurrent form while introducing additional flexibility through gating.
this propagation step, through local nonlinear maps such as multilayer perceptrons (MLPs) or gated linear units (GLUs) [41], which can be viewed as input-dependent gates that modulate the transmitted signal. By contrast, in minGRU, the update gate zt = σ(Wz xt−1 + bz )
makes the recurrence input-dependent, so the nonlinearity appears directly inside the state update, ht = zt ⊙ ht−1 + (1 − zt ) ⊙ (Wh xt−1 + bh ).
Nonlinearity in recurrent architectures
Having introduced SSMs, LRUs, and minGRUs, it is natural to ask where the nonlinearity enters different classes of parallelizable recurrent architectures. Many such models can be written schematically as ht = At ht−1 + Φt (xt−1 , ht−1 ),
(12)
This viewpoint also helps situate other recent architectures. Simplified gated models such as minLSTM [31] follow a similar philosophy to minGRU, where inputdependent gates modulate state propagation. Additionally, selective SSMs, such as Mamba [30], likewise introduce input-dependent updates within the state space framework. By contrast, other approaches dispense with a recurrent hidden state altogether and instead place nonlinearity in an explicit feature map constructed from delayed inputs. Relevant examples include the nonlinear vector autoregression (NVAR) model [42] and next-generation reservoir computing (NGRC) [43]. Concretely, one forms a finite delay vector using k delayed inputs separated by a delay spacing s: xlin (t) = [xt , xt−s , xt−2s , . . . , xt−(k−1)s ],
(13)
and then augments it with polynomial features, for example (p)
xnonlin (t) = ϕp (xlin (t)), ϕp (x) = xi1 xi2 · · · xip 1≤i1 ≤···≤ip ≤k .
(14)
The full feature vector is then ⊤ (2) (p) ot = 1, (xlin (t))⊤ , xnonlin (t)⊤ , . . . , xnonlin (t)⊤ , (15) and the prediction is obtained from a linear readout, e.g. xt+1 = Cot
(16)
or, in the increment form often used in NGRC, xt+1 = xt + Cot .
D.
(11)
(17)
Thus, NVAR/NGRC replaces nonlinear recurrent state evolution by explicit nonlinear feature construction. Overall, a useful categorization principle is not merely whether a model is recurrent, but whether nonlinearity is placed outside the recurrence, inside the recurrent update, or in an explicit feature construction.
(10)
where At governs state propagation and Φt is a function of the previous input and previous hidden state, which injects new information. From this perspective, the LRU keeps the recurrence itself linear, with a fixed diagonal propagator and linear input injection. Nonlinearity is added only outside
E.
Natural encoding of relative distances
Beyond the question of where the nonlinearity is placed, recurrent architectures also possess a natural inductive bias for encoding relative distances [35]. To illustrate this point, it is useful to consider the linearized
4 recursion as a simple example. In this case, the hidden state evolves according to ht = Aht−1 + Bxt ,
(18)
with output yt = Cht . Unrolling the recursion gives ht =
t X
At−j Bxj ,
(19)
CAt−j Bxj .
(20)
j=1
and therefore yt =
t X j=1
The contribution of xj to yt thus depends on the lag t − j through the factor At−j . Equivalently, yt =
t−1 X
CAr Bxt−r ,
(21)
r=0
so that the effective kernel K(r) = CAr B
(22)
depends only on the relative separation r [35]. This shows that, even in this simple linearized setting, the recurrence naturally yields a signal that depends on the relative distances along the sequence in RNNs. This natural relative distance encoding mechanism in RNNs differs from approaches in which geometric information is introduced explicitly through a distancedependent bias. For example, Ref. [20] incorporates Euclidean lattice distances as extra information into the attention weights. By contrast, in a recurrent architecture, including our PSR-NQS framework, distance dependence emerges naturally from the sequential structure itself. F.
Two-dimensional Recurrences
To model two-dimensional systems, we follow Ref. [8] and traverse the Lx × Ly lattice along a fixed zigzag ordering, which preserves autoregressive causality. For a site (i, j), the hidden state receives information from two causal predecessors: one along the horizontal direction and one along the vertical direction. Using the zigzag sampling path convention [8], these predecessor sites are (i − (−1)j , j) and (i, j − 1), with any quantities outside the lattice defined to be zero. (ℓ) Let xi,j denote the input representation at site (i, j)
A generic stacked 2D recurrent layer then updates the state according to (ℓ+1) (ℓ) (ℓ) (ℓ) (ℓ) xi,j , hi,j = F (ℓ) ui,j , hH , hV , (24) where F (ℓ) denotes a generic two-dimensional recurrent cell at depth ℓ. This notation is used to describe the causal 2D architecture independently of the specific re(ℓ) current update. The input ui,j is defined by ( (0) (0) xi−(−1)j ,j , xi,j−1 , ℓ = 0, (ℓ) ui,j = (25) (ℓ) xi,j , ℓ ≥ 1. Thus, at the first layer, the cell takes as input the two causal neighboring physical inputs, while at deeper layers it acts only on the representation at the current site produced by the previous layer. For the two-dimensional benchmarks, we instantiate F (ℓ) using a two-dimensional adaptation of the minGRU introduced in Sec. II C, dubbed 2D minGRU. Here, for a given site (i, j), the horizontal predecessor hidden state (ℓ) hH plays the role of the recurrent state carried along the current row, while the vertical predecessor hidden state (ℓ) hV enters the candidate update as additional causal context. Concretely, the update can be written schematically as (ℓ) (ℓ) zi,j = σ Wz(ℓ) ui,j + b(ℓ) , (26) z (ℓ) (ℓ) (ℓ) (ℓ) (ℓ) e hi,j = Wh ui,j , hV + bh , (27) (ℓ) (ℓ) (ℓ) (ℓ) (ℓ) hi,j . (28) hi,j = zi,j ⊙ hH + 1 − zi,j ⊙ e (ℓ+1)
Finally, the recurrent cell computes its output x̃i,j (ℓ+1)
x̃i,j
(ℓ)
(ℓ)
= C (ℓ) hi,j + D(ℓ) ui,j ,
as
(29)
where the second term denotes a learned linear projection of the inputs. Although the full recurrence is not fully parallelizable through a single parallel scan operation [38], it can still be evaluated efficiently in a row-by-row manner. For fixed j, the recurrence along the snake direction can be computed with a parallel scan over i, reducing the sequential depth within each row from O(L) to O(log L) assuming Lx = Ly = L in the case of a square lattice. The propagation between rows remains sequential, so the recurrent sequential depth is reduced from O(L2 ) to O(L log L), while the total arithmetic work remains proportional to the number of lattice sites, O(L2 ).
(ℓ)
and layer ℓ, and let hi,j be the corresponding hidden state. In this work, we use patches of 2 × 2 spins as (0) initial inputs xi,j to improve parallelization and expressivity [17, 22]. We further define the causal predecessor hidden states as (ℓ)
(ℓ)
hH = hi−(−1)j ,j ,
(ℓ)
(ℓ)
hV = hi,j−1 .
(23)
G.
Residual connections
To facilitate information flow across depth [29, 44], in PSR-NQS, we include a residual pathway between recurrent layers. For the one-dimensional architecture, the input to layer ℓ, called x(ℓ) , is passed to the recurrent layer,
5 whose output we denote by y (ℓ) . The residual contribution is added directly following the GLU [41] network for all layers except the initial layer: ( 0, ℓ < 1, (ℓ+1) (ℓ) (30) xi = GLU(yi ) + (ℓ) xi , ℓ ≥ 1. For the two-dimensional architecture, the residual input to layer ℓ at site (i, j) is defined as ( 0, ℓ ≤ 1, (ℓ) si,j = (31) (ℓ−1) hi,j , ℓ ≥ 2. Then, the residual contribution is added before applying the nonlinear activation [45]: (ℓ+1) (ℓ+1) (ℓ) xi,j = GELU x̃i,j + si,j , (32) (ℓ)
with si,j = 0 for the first two layers. Thus, deeper layers receive both the newly computed recurrent representation and, for ℓ ≥ 2, the hidden state from the previous recurrent layer at the same lattice site. This yields a deep 2D recurrent architecture that preserves autoregressive causality while allowing hierarchical feature propagation across the lattice. III. A.
RESULTS
Runtime benchmarks
To evaluate the practical benefit of the parallel scan implementation, we compare the runtime of the sequential recurrence (sequential) and parallel scan recurrence (parallel) modes for two representative architectures used in this work: the 1D LRU and the 2D minGRU. The goal of these benchmarks is not to compare different models against one another, but rather to quantify the speedup obtained by replacing the sequential recurrent evaluation with its parallel counterpart whenever available. For each architecture, we measure the wall-clock time of a single training step in both modes under the same hardware and software conditions. Each training step consists of autoregressive sampling of computational basis configurations, calculating the local energy, computing the gradient, and updating the parameters. In order to make the comparison as direct as possible, the model size, batch size, and input dimensions are kept fixed between the two modes. In the 1D LRU model, the recurrence can be evaluated either sequentially or by means of parallel scan, allowing for a direct comparison between the two execution modes. In the 2D setting, although the full recurrence is not globally parallelizable, the 2D minGRU can still be evaluated in a row-by-row manner, where the recurrence along each row is carried out using parallel scan, while the recurrence across different rows remains sequential.
FIG. 1. Runtime per training step comparison between the parallel scan and sequential versions of (a) the 1D LRU for the 1D TFIM with a system size N = L and (b) the 2D minGRU with 2 × 2 patchesfor the square lattice Heisenberg antiferromagnet with a system size N = L × L. Dotted lines showing the ideal scaling laws are included as guides to the eye. We observe that the parallel scan implementation is faster than the sequential recurrent mode. The runtime comparisons were conducted on one H100 GPU.
This provides a practical acceleration while preserving the autoregressive structure. The benchmark results reported in Fig. 1 show that the parallel implementation yields a clear reduction in training step runtime for both the 1D LRU and the 2D minGRU. Our asymptotic scaling is consistent with the predicted scaling O(L log L) in 1D and O(L3 log L) in 2D. Overall, our results show that the PSR-NQS formulation leads to a tangible practical speedup in a VMC setting. We note, however, that the attainable parallel scan speedup is ultimately limited by GPU memory, since the parallel formulation trades reduced sequential depth for increased memory traffic, and requires being able to store the entire sequence of hidden states [28, 29].
6 B.
One-dimensional Transverse-field Ferromagnetic Ising model
As a first benchmark, we consider the one-dimensional TFIM, with Hamiltonian Ĥ = −
N −1 X i=1
z σiz σi+1 −h
N X
σix ,
(33)
i=1
where h is the transverse-field strength and N is the number of spins. Here, σix,z are Pauli operators acting on site i. This model provides a standard testbed for a variational wave function ansatz, as it is exactly solvable and exhibits a competition between ferromagnetic ordering and quantum fluctuations induced by the transverse field h, with a critical point at h = 1 [46, 47], which we focus on in this study. The variational state is optimized within the VMC framework [1] by minimizing the expectation value of the Hamiltonian Ĥ, E(θ) = ⟨Ψθ |Ĥ|Ψθ ⟩,
(34)
within our normalized PSR-NQS |Ψθ ⟩. In practice, expectation values are estimated from samples σ drawn exactly and autoregressively from the probability distribution |Ψθ (σ)|2 , where the amplitude Ψθ (σ) is computed following the formalism in Ref. [8]. The variational parameters are updated iteratively with Adam optimizer using stochastic estimates of the energy gradient [8, 48]. This procedure allows our ansatz to approximate the ground state of our Hamiltonians of interest through direct energy minimization. We train a deep LRU with three layers, with each block consisting of a 1D LRU (with dh = 64), followed by a single-layer GLU feed-forward network [41] acting independently on each site. The ansatz is trained using the iterative retraining technique [9, 15, 49], in which the model is first trained at a small system size, and then trained further at progressively larger system sizes. We then perform a finite-size extrapolation to the groundstate energy per site in the thermodynamic limit by fitting the variational energy densities obtained at each system size to the scaling form: eNQS (N ) = e∞ NQS +
The runtime increases with system size, from about 8 ms per step at N = 6 to about 0.9 s per step at N = 256. Overall, the full set of calculations required slightly less than nine hours on a single H100 GPU, showing that the iterative retraining strategy yields high-accuracy energies at moderate computational cost. More details about our hyperparameters are provided in Appendix A.
a2 a3 a1 + 2 + 3. N N N
(35)
We report details of the numerical fit in Appendix B. Tab. I demonstrates that the 1D LRU ansatz accurately reproduces the exact open-boundary TFIM energies at the critical point h = 1. The relative errors remain below 10−4 for all system sizes up to N = 256, with several cases reaching the 10−6 level or better. Furthermore, the energy density estimate in the thermodynamic limit is within a small relative error of 3.11 × 10−5 compared to the exact value. The non-monotonic dependence of the error on N suggests that the remaining discrepancies are dominated by stochastic optimization noise rather than by a systematic loss of expressivity.
C.
Square Lattice Heisenberg model
We benchmark our variational ansatz on the spin- 12 two-dimensional antiferromagnetic Heisenberg model on an L × L square lattice with open boundary conditions (OBC). The Hamiltonian reads X Ĥ = Si · Sj , (36) ⟨i,j⟩
where ⟨i, j⟩ runs over nearest-neighbor bonds of the square lattice, and Si = (Six , Siy , Siz ) denotes a spin- 12 operator at site i with components Siα = 12 σiα (α = x, y, z) in terms of the Pauli matrices σiα . This model is particularly useful as a benchmark for variational wave functions for two main reasons. First, it has long served as a canonical testbed for ground-state methods in strongly correlated quantum matter [15]. In fact, it is an interacting two-dimensional quantum system with gapless Goldstone modes and long-ranged spin correlations associated with Néel order, which places strong demands on expressivity beyond short-range physics [53, 54]. Second, because the model becomes sign-free on the bipartite square lattice, after applying a Marshall sign [55, 56], high-precision stochastic series expansion quantum Monte Carlo (QMC) estimates are available and serve as a reliable reference for assessing the absolute accuracy of variational energies [15, 52]. Tab. II compares our results to representative variational baselines: projected-pair entangled states (PEPS) [51], PixelCNN [10], and a 2D tensorized RNN (2D TRNN) [15], in addition to QMC reference values [52]. Our 2D minGRU architecture consists of 6 stacked layers, where each layer has a hidden dimension dh = 512. To further improve our accuracy, we impose the c4v point group symmetry [8, 57]. On the 10 × 10 system, our 2D minGRU achieves a competitive energy with the best reported RNN variational energy and lies within 1.9 × 10−5 per site of the QMC reference. On the larger 16 × 16 cluster, our 2D minGRU yields the lowest variational energy in our comparison, improving upon PEPS and PixelCNN and approaching the QMC reference within 2.8 × 10−5 per site. Overall, the 2D minGRU remains accurate despite using a linearized recursion and forgoing tensorization, used in 2D TRNNs to increase expressivity [15, 57–59], while simultaneously improving computational efficiency through faster recurrent updates. We now turn our attention to the performance of the iteratively retrained 2D minGRU, with 3 layers and a hid-
7 TABLE I. A comparison between 1D LRU energies per site e1DLRU (N ) (with the one-standard error uncertainties on the last digits in parentheses) and the exact 1D TFIM per-site-energies for OBC at h = 1, following the formula eexact (N ) = [1 − csc(π/(2(2N + 1)))]/N [50]. The relative error with respect to the exact reference energies is reported in units of 10−5 . We also a3 a2 report an extrapolation to the thermodynamic limit by fitting the variational energies to the model e(N ) = e∞ + aN1 + N 2 + N3 . The cumulative training time, hereafter, is reported in hh:mm:ss format. We also report our runtimes for training the 1D LRU using a single H100 95GB GPU. VMC is performed using 1024 samples at each iteration, and the final energies are computed at the end of each size’s training procedure using 1,024,000 samples. For reference, this model consists of 3 layers with d = dh = 64 and has around 75,000 variational parameters. N
e1DLRU (N )
eexact (N ) Rel. err. (×10−5 )
6 −1.2160344(3) −1.2160383 8 −1.2297437(5) −1.2297439 10 −1.2381549(7) −1.2381490 12 −1.2438271(9) −1.2438309 16 −1.251015(1) −1.2510242 24 −1.258337(1) −1.2583213 32 −1.262069(2) −1.2620098 48 −1.265666(2) −1.2657254 64 −1.267473(2) −1.2675934 96 −1.269452(1) −1.2694685 128 −1.270408(1) −1.2704086 192 −1.271305(1) −1.2713505 256 −1.271797(1) −1.2718221 ∞ −1.2731999(8) −1.2732395
Time per training Number of Cumulative step (ms) training steps training time
0.33 0.01 0.47 0.31 0.72 1.28 4.69 4.70 9.52 1.33 0.03 3.58 2.00 3.11
8 9 10 12 14 20 28 49 71 150 260 520 900 —
55000 39000 30000 24000 18000 15000 15000 15000 15000 15000 15000 15000 15000 —
00:07:20 00:13:15 00:18:23 00:23:02 00:27:19 00:32:26 00:39:19 00:51:34 01:09:15 01:47:21 02:51:13 05:01:15 08:45:18 —
TABLE II. Ground-state energy per site E/N estimates for the 2D Heisenberg model on the square lattice with OBC. The best energies among the variational energies (columns second through fifth from the left) are bolded. QMC energies are provided as a reference. Numbers in parentheses denote one-standard-error uncertainties in the last digits. Note that our 10×10 calculation took around 4.5 GPU days, whereas 16 × 16 calculation was conducted in a time frame of 29 GPU days using a single L40S 48GB GPU. For reference, this 2D minGRU model has 6 layers and a hidden dimension dh = 512 with around 7.1 million variational parameters. Lattice
PEPS [51]
PixelCNN [10] 2D TRNN [15] 2D minGRU (ours)
10 × 10 −0.628601(2) −0.628627(1) −0.628656(9) 16 × 16 −0.643391(3) −0.643448(1) —
den dimension size dh = 256, as summarized in Fig. 2. Starting from initial training on a 6×6 lattice, we reached larger sizes, up to L × L = 52 × 52, through iterative retraining, initializing each new system size from the optimized parameters of the previous size. In this way, a substantial part of the optimization is effectively carried out on smaller lattices, which provide favorable initial conditions for training at larger sizes and a natural setting for conducting finite-size scaling studies without restarting from scratch at each system size [9, 15, 22, 49, 57, 60]. The results show consistently good accuracy over a wide range of lattice sizes, with relative errors of order 10−4 compared to QMC, while all calculations were performed on a single A100 GPU. More details about our simulation hyperparameters are provided in Appendix A. To obtain an estimate of the thermodynamic limit ground-state energy per site, we perform a finite-size extrapolation of the 2D minGRU energies using the following fitting model [52]: eNQS (L) = e∞ NQS +
a1 a2 a3 + 2 + 3. L L L
(37)
−0.628637(4) −0.643504(3)
QMC [52] −0.6286561(2) −0.6435317(2)
Our thermodynamic-limit estimate e∞ NQS , reported in the inset of Fig. 2, is consistent with the corresponding QMC extrapolation within a relative error of 7.5 × 10−4 . The full numerical data used in this finite-size analysis are provided in Appendix C. These results indicate that the recurrent ansatz remains competitive even in a large-scale regime that is already highly challenging for neural quantum states without deploying significant GPU resources. In particular, the 52 × 52 calculation is, to our knowledge, among the largest lattice sizes reported for NQS using a single A100 GPU, improving on the recent 42 × 42 benchmark obtained with a vision transformer (ViT)-based approach [20] and the 40 × 40 benchmark reported for patched transformers [17]. As a last note, we would like to highlight that recent transformer NQS constructions [20] introduce an explicit decaying-correlation inductive bias to stabilize large-system optimization. Our approach does not require imposing such a bias by hand. Here distancedependent correlations arise naturally from the recurrent
8 2D minGRU QMC Interpolated QMC
−0.61
−0.63 8
−0.64
6
[10−4 ]
E/N
−0.62
−0.65
2
−0.66 −0.67
4
0 0.00
0.05
0.10
0.15
1/L
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
0.16
1/L FIG. 2. Finite-size scaling of the ground-state energy per spin E/N for the square lattice Heisenberg model with OBC for L = 6 to 52, corresponding to N = L2 spins. The 2D minGRU results are shown as circles. QMC reference values listed in Tab. VI are shown as filled triangles, while additional interpolated QMC values are shown as open triangles. The points at 1/L = 0 denote the corresponding infinitesize extrapolations, while the dotted lines serve as guides to the eye. The inset shows the relative error, defined as ϵ = [E2DminGRU − EQMC ]/|EQMC |, plotted in units of 10−4 . The relative error remains small for all system sizes, staying at the level of a few 10−4 , and remains below 10−3 even for the infinite-size extrapolation. Error bars are shown, but are smaller than the marker sizes.
architecture itself as highlighted in Ref. [35]. IV.
CONCLUSION
In this work, we have shown that PSR-NQS provide a natural and efficient framework for scalable simulations of many-body systems. A central advantage of the recurrent construction is its inductive bias toward relative distance, without requiring this structure to be imposed explicitly, as is commonly done in transformer-based architectures. This makes recurrent models particularly appealing for problems where distance and correlation structure play an essential role. From a computational perspective, the resulting VMC optimization of PSR-NQS is favorable. In particular, the training time scales as O(L log L) in one dimension and as O(L3 log L) in two dimensions, allowing us to reach larger system sizes more efficiently than regular RNN wave functions with O(L2 ) and O(L4 ) scaling in oneand two-spatial dimensions, respectively. Numerically, our two-dimensional minGRU ansatz yields state-of-theart ground-state energies for the square-lattice Heisenberg model on benchmark systems such as the 10 × 10 and 16 × 16 lattices. Moreover, by combining the ansatz with iterative retraining, we are able to push VMC simulations up to 52 × 52 sites, while maintaining strong
agreement with available QMC results, all on a single GPU. Overall, an important practical message of this work is that large-scale NQS simulations do not necessarily require heavy computational infrastructure. The total computational budget used for training in all our simulations is estimated to be around 1400 GPU hours, which remains modest in comparison to recent large-scale transformer studies [20]. Although the present framework could be pushed further through newer hardware and more aggressive parallel GPU implementations, this is not essential to obtain accurate results. Rather, our results show that accurate and scalable neural quantum state calculations can already be achieved with relatively light computational resources, making such simulations accessible to a broader range of research groups. Extending these conclusions to more general settings, such as fermionic systems, remains an important direction for future work. ACKNOWLEDGMENTS
We are grateful to Schuyler Moss for helpful discussions. E.M., E.K., and R.S. were supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award Number DESC0022311. M.H. and M.K. acknowledge support from the Natural Sciences and Engineering Research Council of Canada (NSERC). Some of the ideas for this work were formed during our time at the Kavli Institute for Theoretical Physics (KITP), which is supported in part by the NSF grant PHY-2309135 and by the HeisingSimons Foundation. Computer simulations were made possible thanks to University of Waterloo’s Math Faculty Computing Facility (MFCF) and San José State University’s GPU cluster acquired through the National Science Foundation Grant No. OAC-2430291. While writing this manuscript, we became aware of another paper that also used parallel scan techniques through DysonNet [61]. We arrived at our results independently and are excited to see further applications of parallel scan to a non-autoregressive architecture in the NQS literature. We would like to emphasize that the focus of our paper is on optimized autoregressive wave functions using parallel scan techniques, in both one and two spatial dimensions. CODE AVAILABILITY
Our implementation of the presented methods and all scripts needed to reproduce our results in this manuscript are openly available on GitHub https://github.com/ ParallelScan-RNNs/PSR-NQS.
9 exact scaling parameters:
Appendix A: Hyperparameters
In this appendix, we summarize all training settings for reproducibility purposes. The one-dimensional TFIM simulation hyperparameters are provided in Tab. III. For two-dimensional simulations, cold-start runs for N = 102 , 162 are listed in Tab. IV. Iterative retraining is also performed starting from a cold-start simulation at L = 6, with settings summarized in Tab. V. For learning rate decay, we use −1
η(t) = η0 (1 + t/δ0 )
,
(A1)
and for iterative retraining at L > L0 we set the number of training steps to [15, 49] l m Nsteps (L) = s Ce−r(L−L0 ) + F . (A2) Here, L is the linear lattice size, s is an overall scale factor, C and r set the amplitude and decay rate of the size-dependent term, L0 is a reference size, and F is the large-L offset. The stage protocol inspired by Refs. [15, 49, 57] is as follows. At the initial size L = L0 = 6, training is split into three consecutive stages: (i) stage 1 uses no symmetries (nosym) with fixed learning rate ηstage1 ; (ii) stage 2 keeps nosym but switches to the decaying schedule in Eq. (A1) for an additional number of steps denoted as Nstage2 . (iii) Stage 3 switches to c4v symmetry and continues with Eq. (A1) for Nstage3 . For all larger sizes L > L0 , training uses a single c4v symmetry training stage with fixed learning rate ηiter , and the training budget is set by Eq. (A2). Here ηiter denotes the fixed learning rate used for iterative retraining at L > L0 , while ηstage1 is the fixed learning rate used in the initial stage at L = L0 . The parameters η0,stage and δ0 define the decaying schedule in Eq. (A1) for later stages at L = L0 .
4 ≈ −1.27323954, π π−2 aexact = ≈ 0.36338023, 1 π π aexact =− ≈ −0.13089969, 2 24 π aexact = ≈ 0.06544985. 3 48 eexact =− ∞
We remark that although the finite-size fit gives a remarkably good estimate of the per-site energy in the thermodynamic limit, e∞ , with relative error on the order of 10−5 , the other coefficients are much more difficult to fit correctly with a truncated series, due to the illconditioned nature of the fitting problem, which involves powers of 1/N . Appendix C: Finite-size scaling on the square lattice Heisenberg model
In this appendix, we provide the numerical data used for the finite-size scaling analysis shown in Fig. 2. The results include the 2D minGRU variational energies, the corresponding QMC reference values [52], the relative error, and the training cost for each lattice size. The relative error is defined as ϵ=
ENQS − EQMC . |EQMC |
The full data are reported in Tab. VI. The points labeled L = ∞ correspond to independent infinite-size extrapolations of the finite-size 2D minGRU (NQS) and QMC data. Both the 2D minGRU and QMC finite-size data are fit to the open-boundary scaling form e(L) = e∞ +
Appendix B: Finite-size scaling on the 1D Transverse-Field Ising model
a2 a3 a1 + 2 + 3, L L L
(C1)
where e(L) = E(L)/N . The fitted parameters are We provide numerical data for the thermodynamic limit extrapolation presented in Tab. I. The variational energy densities are fit to the scaling form: e(N ) = e∞ +
a2 a3 a1 + 2 + 3, N N N
(B1)
where e(N ) = E(N )/N . We obtain the fitted parameters e∞ = −1.2731999(8), a2 = −0.1179(7),
a1 = 0.36207(5), a3 = 0.021(3).
The exact form of the energy density for the 1D TFIM with open boundaries at the critical point is given by eexact (N ) = [1 − csc(π/(2(2N + 1)))]/N [50], and can be expanded about N = ∞ via a Laurent series, giving the
eQMC = −0.6694585(2), ∞
eNQS = −0.668938(6), ∞
aQMC = 0.42583(1), 1
aNQS = 0.4150(3), 1
= −0.1734(2), aQMC 2
aNQS = −0.090(4), 2
aQMC = −0.0473(6), 3
aNQS = −0.27(2). 3
The 2D minGRU and QMC extrapolations show close agreement over the full finite-size range considered. In particular, the fitted thermodynamic limit energies differ by about 5 × 10−4 per spin, corresponding to a relative difference below 10−3 . For the finite-size comparison in Fig. 2, however, we use QMC data with the same open boundary conditions as the 2D minGRU calculations and fit them using the OBC scaling form in Eq. (C1). This OBC fit is used to
10 TABLE III. Hyperparameters for 1D LRU runs with iterative retraining from a cold start at N = 6. Simulation
Hyperparameter
Value
Architecture 1D LRU with skip-connections Number of layers 3 Hidden/model dimensions dh = d = 64 Number of samples 1024 N = 6 to 256 Learning rate η = 1 × 10−4 Adam optimizer parameters β1 = 0.9, β2 = 0.999, ϵ = 10−8 Eq. (A2) constants s = 1.0, r = 0.25, L0 = 6, C = 40000, F = 15000 (Nsteps is then rounded to the nearest multiple of 1000)
TABLE IV. Hyperparameters for 2D minGRU runs with cold starts. Simulation Hyperparameter
Value
N = 102
Architecture 2D patched minGRU with skip-connections Number of layers 6 Hidden dimension size dh = 512 Patch size (px , py ) = (2, 2) Symmetry c4v Number of samples 200 Training iterations 150000 Learning rate η(t) = 5 × 10−4 (1 + t/5000)−1
N = 162
Architecture 2D patched minGRU with skip-connections Number of layers 6 Hidden dimension size dh = 512 Patch size (px , py ) = (2, 2) Symmetry c4v Number of samples 200 Training iterations 150000 Learning rate η(t) = 5 × 10−4 (1 + t/5000)−1
compare the finite-size trends and to define the relative errors shown in the inset of Fig. 2. However, we ac-
knowledge that high-precision QMC extrapolation with periodic boundary conditions provides a more reliable estimate [52].
[1] F. Becca and S. Sorella, Quantum Monte Carlo Approaches for Correlated Systems (Cambridge University Press, 2017). [2] G. Carleo and M. Troyer, Solving the quantum manybody problem with artificial neural networks, Science 355, 602 (2017). [3] H. Lange, A. Van de Walle, A. Abedinnia, and A. Bohrdt, From architectures to applications: a review of neural quantum states, Quantum Science and Technology 9, 040501 (2024). [4] M. Medvidović and J. R. Moreno, Neural-network quantum states for many-body physics, The European Physical Journal Plus 139, 631 (2024). [5] A. Dawid, J. Arnold, B. Requena, A. Gresch, M. Plodzień, K. Donatella, K. A. Nicoli, P. Stornati, R. Koch, M. Büttner, R. Okula, G. Muñoz-Gil, R. A. Vargas-Hernández, A. Cervera-Lierta, J. Carrasquilla, V. Dunjko, M. Gabrié, P. Huembeli, E. van 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). [6] K. Choo, T. Neupert, and G. Carleo, Two-dimensional frustrated J1 -J2 model studied with neural network quantum states, Physical Review B 100, 125124 (2019). [7] D. Luo and B. K. Clark, Backflow transformations via neural networks for quantum many-body wave functions, Physical Review Letters 122, 10.1103/physrevlett.122.226401 (2019). [8] M. Hibat-Allah, M. Ganahl, L. E. Hayward, R. G. Melko, and J. Carrasquilla, Recurrent neural network wave functions, Physical Review Research 2, 10.1103/physrevresearch.2.023358 (2020). [9] C. Roth, Iterative retraining of quantum spin models using recurrent neural networks (2020), arXiv:2003.06228 [physics.comp-ph]. [10] O. Sharir, Y. Levine, N. Wies, G. Carleo, and A. Shashua, Deep autoregressive models for the efficient variational simulation of many-body quantum systems, Phys. Rev. Lett. 124, 020503 (2020). [11] C. Roth, A. Szabó, and A. MacDonald, High-accuracy variational Monte Carlo for frustrated magnets with deep
11 TABLE V. Hyperparameters for 2D minGRU iterative retraining initialized from a cold-start at N = 62 . System size
Hyperparameter
Value
N = 62 to 522 (all iterative runs)
Architecture 2D minGRU Number of layers 3 Hidden/model size dh = d = 256 Patch size (px , py ) = (2, 2) Number of samples (train) 200 Stage-1 learning rate ηstage1 = 5 × 10−4 Stage learning-rate base η0,stage = 5 × 10−4 Decay-scale parameter δ0 = 5000 Stage 1 schedule params Nstage1 = 201000 Stage 2 iterations Nstage2 = 76000 Stage 3 iterations Nstage3 = 101000 Eq. (A2) constants r = 0.25, L0 = 6, C = 101000 Final evaluation symmetry c4v , samples = 100000 Symmetry Large L learning rate Eq. (A2) exponent Eq. (A2) offset
c4v ηiter = 5 × 10−5 s = 4.0 F = 2000
Symmetry Large L learning rate N = 18 , 20 , 22 , 24 , 26 , 28 , 30 , 32 Eq. (A2) exponent Eq. (A2) offset
c4v ηiter = 1 × 10−4 s = 1.0 F = 2000
Symmetry Large L learning rate Eq. (A2) exponent Eq. (A2) offset
c4v ηiter = 1 × 10−4 s = 1.0 F = 1000
Symmetry Large L learning rate Eq. (A2) exponent Eq. (A2) offset
c4v ηiter = 1 × 10−4 s = 1.0 F = 500
2
2
2
2
2
N = 6 , 8 , 10 , 12 , 14 , 16
2
2
2
2
2
2
N = 36 , 40 , 44
2
N = 462 , 482 , 502 , 522
2
2
2
2
2
neural networks (2023). [12] Y.-H. Zhang and M. Di Ventra, Transformer quantum state: A multipurpose model for quantum manybody problems, Physical Review B 107, 10.1103/physrevb.107.075147 (2023). [13] T. D. Barrett, A. Malyshev, and A. I. Lvovsky, Autoregressive neural-network wavefunctions for ab initio quantum chemistry, Nature Machine Intelligence 4, 351 (2022). [14] A. Chen and M. Heyl, Empowering deep neural quantum states through efficient optimization, Nature Physics 20, 1476–1481 (2024). [15] M. S. Moss, R. Wiersema, M. Hibat-Allah, J. Carrasquilla, and R. G. Melko, Leveraging recurrence in neural network wavefunctions for large-scale simulations of Heisenberg antiferromagnets on the square lattice, Phys. Rev. B 112, 134450 (2025). [16] R. Rende, L. L. Viteritti, L. Bardone, F. Becca, and S. Goldt, A simple linear algebra identity to optimize large-scale neural network quantum states, Communications Physics 7, 10.1038/s42005-024-01732-4 (2024). [17] K. Sprague and S. Czischek, Variational monte carlo with large patched transformers, Communications Physics 7, 10.1038/s42005-024-01584-y (2024). [18] D. S. Kufel, J. Kemp, D. Vu, S. M. Linsel, C. R. Laumann, and N. Y. Yao, Approximately symmetric neural
networks for quantum spin liquids, Physical Review Letters 135, 10.1103/pgnx-11ph (2025). [19] Y. Nomura and M. Imada, Dirac-type nodal spin liquid revealed by refined quantum many-body solver using neural-network wave function, correlation ratio, and level spectroscopy, Phys. Rev. X 11, 031034 (2021). [20] L. L. Viteritti, R. Rende, S. Sachdev, and G. Carleo, Approaching the thermodynamic limit with neural-network quantum states (2026), arXiv:2602.02665 [cond-mat.strel]. [21] Y. Gu, W. Li, H. Lin, B. Zhan, R. Li, Y. Huang, D. He, Y. Wu, T. Xiang, M. Qin, L. Wang, and D. Lv, Solving the Hubbard model with neural quantum states (2025), arXiv:2507.02644 [cond-mat.str-el]. [22] M. Hibat-Allah, E. Merali, G. Torlai, R. G. Melko, and J. Carrasquilla, Recurrent neural network wave functions for Rydberg atom arrays on kagome lattice, Communications Physics 8, 10.1038/s42005-025-02226-7 (2025). [23] H. Lange, A. Chen, A. Georges, F. Grusdt, A. Bohrdt, and C. Roth, Simulating superconductivity in mixeddimensional t∥ -J∥ -J⊥ bilayers with neural quantum states (2026), arXiv:2602.10091 [cond-mat.str-el]. [24] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, Attention is all you need, in Advances in Neural Information Processing Systems, Vol. 30, edited by I. Guyon,
12 TABLE VI. A comparison between 2D minGRU energies per site (with 3 layers and c4v symmetry trained using the iterative retraining technique) and QMC (SSE) data [52] on the 2D Heisenberg model with OBC. The relative error ϵ, with respect to a2 a3 QMC, is reported in units of 10−4 . We perform a fitting using the model e(L) = e∞ + aL1 + L 2 + L3 to obtain the thermodynamic limit energy extrapolation of the 2D minGRU and QMC energies. We also report our runtimes for training the 2D minGRU using a single A100 80GB GPU. For reference, this PSR-NQS model has around 800,000 variational parameters. N
2D minGRU
QMC
62 −0.603519(2) −0.6035222(2) 82 −0.619025(5) −0.6190371(2) 102 −0.628605(6) −0.6286561(2) 122 −0.635127(7) −0.6352007(2) 142 −0.639842(6) −0.6399410(2) 162 −0.643396(6) −0.6435317(2) 182 −0.646186(6) −0.6463451(4) 202 −0.648401(5) −0.6486091(4) 222 −0.650253(5) −0.6504689(4) 242 −0.651812(5) −0.6520251(4) 262 −0.653129(5) −0.6533456(4) 282 −0.654261(4) −0.6544806(4) 302 −0.655241(4) 322 −0.656089(4) −0.6563289(4) 362 −0.657523(4) 402 −0.658657(4) 442 −0.659565(4) 462 −0.659941(4) 482 −0.660323(4) −0.6606690(4) 502 −0.660649(3) 522 −0.660967(4) ∞ −0.668938(6) −0.6694585(2)
Rel. err. (×10−4 ) Time per training step (s) Training steps Cumulative training time 0.1 0.2 0.8 1.2 1.5 2.1 2.5 3.2 3.3 3.3 3.3 3.4 3.7 5.2 7.8
U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017). [25] E. Martin and C. Cundy, Parallelizing linear recurrent neural nets over sequence length (2018), arXiv:1709.04057 [cs.NE]. [26] A. Gu, T. Dao, S. Ermon, A. Rudra, and C. Re, Hippo: Recurrent memory with optimal polynomial projections (2020), arXiv:2008.07669 [cs.LG]. [27] A. Gu, I. Johnson, K. Goel, K. Saab, T. Dao, A. Rudra, and C. Ré, Combining recurrent, convolutional, and continuous-time models with linear state-space layers (2021), arXiv:2110.13985 [cs.LG]. [28] A. Gu, K. Goel, and C. Ré, Efficiently modeling long sequences with structured state spaces (2022), arXiv:2111.00396 [cs.LG]. [29] A. Orvieto, S. L. Smith, A. Gu, A. Fernando, C. Gulcehre, R. Pascanu, and S. De, Resurrecting recurrent neural networks for long sequences (2023), arXiv:2303.06349 [cs.LG]. [30] A. Gu and T. Dao, Mamba: Linear-time sequence modeling with selective state spaces (2024), arXiv:2312.00752 [cs.LG]. [31] L. Feng, F. Tung, M. O. Ahmed, Y. Bengio, and H. Hajimirsadeghi, Were RNNs all we needed? (2024), arXiv:2410.01201 [cs.LG]. [32] M. Beck, K. Pöppel, M. Spanring, A. Auer, O. Prudnikova, M. Kopp, G. Klambauer, J. Brandstetter, and S. Hochreiter, xLSTM: Extended long short-term memory (2024), arXiv:2405.04517 [cs.LG]. [33] A. Katharopoulos, A. Vyas, N. Pappas, and F. Fleuret,
0.1 0.2 0.5 1.0 1.8 3.3 5.4 8.2 11.6 16.6 22.5 31.6 40.5 57.0 89.3 139.1 196.1 229.3 274.6 331.4 405.6 -
378000 253039 156624 98145 62676 41163 7029 5050 3850 3123 2681 2413 2251 2152 1056 1021 1008 505 503 502 502 -
10:30:00 24:33:28 46:18:40 73:34:25 104:54:42 142:38:40 153:11:16 164:41:26 177:05:46 191:29:48 208:15:10 229:26:01 254:45:27 288:49:51 315:01:32 354:28:33 409:23:01 441:32:58 479:55:02 526:07:45 582:41:16 -
Transformers are rnns: fast autoregressive transformers with linear attention, in Proceedings of the 37th International Conference on Machine Learning, ICML’20 (JMLR.org, 2020). [34] B. Peng, E. Alcaide, Q. Anthony, A. Albalak, S. Arcadinho, S. Biderman, H. Cao, X. Cheng, M. Chung, M. Grella, K. K. GV, X. He, H. Hou, J. Lin, P. Kazienko, J. Kocon, J. Kong, B. Koptyra, H. Lau, K. S. I. Mantri, F. Mom, A. Saito, G. Song, X. Tang, B. Wang, J. S. Wind, S. Wozniak, R. Zhang, Z. Zhang, Q. Zhao, P. Zhou, Q. Zhou, J. Zhu, and R.-J. Zhu, RWKV: Reinventing RNNs for the transformer era (2023), arXiv:2305.13048 [cs.CL]. [35] A. B. Ayub, A. M. Aboussalah, and M. Hibat-Allah, Geometry-induced long-range correlations in recurrent neural network quantum states (2026), arXiv:2604.08661 [quant-ph]. [36] A. Tustin, A method of analysing the behaviour of linear systems in terms of time series, Journal of the Institution of Electrical Engineers - Part IIA: Automatic Regulators and Servo Mechanisms 94, 130 (1947), https://digitallibrary.theiet.org/doi/pdf/10.1049/ji-2a.1947.0020. [37] M. Zhang, K. K. Saab, M. Poli, T. Dao, K. Goel, and C. Ré, Effectively modeling time series with simple discrete state spaces (2023), arXiv:2303.09489 [cs.LG]. [38] G. E. Blelloch, Prefix sums and their applications (1990). [39] M. Schmitt and M. Heyl, Quantum many-body dynamics in two dimensions with artificial neural networks, Phys. Rev. Lett. 125, 100503 (2020). [40] K. Cho, B. van Merrienboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio, Learning
13 phrase representations using RNN encoder-decoder for statistical machine translation (2014), arXiv:1406.1078 [cs.CL]. [41] N. Shazeer, GLU variants improve transformer (2020), arXiv:2002.05202 [cs.LG]. [42] J. H. Adler, S. Hocking, X. Hu, and S. Islam, Physics-informed nonlinear vector autoregressive models for the prediction of dynamical systems (2024), arXiv:2407.18057 [math.DS]. [43] D. J. Gauthier, E. Bollt, A. Griffith, and W. A. S. Barbosa, Next generation reservoir computing, Nature Communications 12, 10.1038/s41467-021-25801-2 (2021). [44] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition (2015), arXiv:1512.03385 [cs.CV]. [45] D. Hendrycks and K. Gimpel, Gaussian error linear units (gelus) (2023), arXiv:1606.08415 [cs.LG]. [46] P. Pfeuty, The one-dimensional ising model with a transverse field, Annals of Physics 57, 79 (1970). [47] G. B. Mbeng, A. Russomanno, and G. E. Santoro, The quantum ising chain for beginners, SciPost Physics Lecture Notes 10.21468/scipostphyslectnotes.82 (2024). [48] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2017), arXiv:1412.6980 [cs.LG]. [49] M. S. Moss, R. Wiersema, M. Hibat-Allah, J. Carrasquilla, and R. G. Melko, Leveraging recurrence in neural network wavefunctions for large-scale simulations of Heisenberg antiferromagnets on the triangular lattice, Phys. Rev. B 112, 134449 (2025). [50] M. Campostrini, A. Pelissetto, and E. Vicari, Quantum ising chains with boundary fields, Journal of Statistical Mechanics: Theory and Experiment 2015, P11015 (2015). [51] W.-Y. Liu, S.-J. Dong, Y.-J. Han, G.-C. Guo, and L. He, Gradient optimization of finite projected entan-
gled pair states, Physical Review B 95, 10.1103/physrevb.95.195154 (2017). [52] A. W. Sandvik, High-precision ground state parameters of the two-dimensional spin-1/2 Heisenberg model on the square lattice (2026), arXiv:2601.20189 [cond-mat.str-el]. [53] P. W. Anderson, An approximate quantum theory of the antiferromagnetic ground state, Phys. Rev. 86, 694 (1952). [54] S. Chakravarty, B. I. Halperin, and D. R. Nelson, Twodimensional quantum Heisenberg antiferromagnet at low temperatures, Phys. Rev. B 39, 2344 (1989). [55] W. Marshall, Antiferromagnetism, Proceedings of the Royal Society of London Series A 232, 48 (1955). [56] M. A. Shamim, M. M. R. Raj, M. Hibat-Allah, and P. T. Araujo, Graph-theoretic analysis of phase optimization complexity in variational wave functions for heisenberg antiferromagnets (2026), arXiv:2602.04943 [condmat.str-el]. [57] M. Hibat-Allah, R. G. Melko, and J. Carrasquilla, Supplementing recurrent neural network wave functions with symmetry and annealing to improve accuracy (2024), arXiv:2207.14314 [cond-mat.dis-nn]. [58] M. Hibat-Allah, E. M. Inack, R. Wiersema, R. G. Melko, and J. Carrasquilla, Variational neural annealing, Nature Machine Intelligence 3, 952–961 (2021). [59] D. Wu, R. Rossi, F. Vicentini, and G. Carleo, From tensor-network quantum states to tensorial recurrent neural networks, Physical Review Research 5, 10.1103/physrevresearch.5.l032001 (2023). [60] D. Luo, Z. Chen, K. Hu, Z. Zhao, V. M. Hur, and B. K. Clark, Gauge-invariant and anyonic-symmetric autoregressive neural network for quantum lattice models, Phys. Rev. Res. 5, 013216 (2023). [61] L. Winter and A. Nunnenkamp, DysonNet: Constanttime local updates for neural quantum states (2026), arXiv:2603.11189 [quant-ph].