Neural-Network Quantum States for Spin-1 Systems: Spin-Basis and Parameterization Effects on Compactness of Representations

Abstract

Neural network quantum states (NQS) have been widely applied to spin-1/2 systems, where they have proven to be highly effective. The application to systems with larger on-site dimension, such as spin-1 or bosonic systems, has been explored less and predominantly using spin-1/2 Restricted Boltzmann Machines (RBMs) with a one-hot/unary encoding. Here, we propose a more direct generalization of RBMs for spin-1 that retains the key properties of the standard spin-1/2 RBM, specifically trivial product states representations, labeling freedom for the visible variables and gauge equivalence to the tensor network formulation. To test this new approach, we present variational Monte Carlo (VMC) calculations for the spin-1 anti-ferromagnetic Heisenberg (AFH) model and benchmark it against the one-hot/unary encoded RBM demonstrating that it achieves the same accuracy with substantially fewer variational parameters. Furthermore, we investigate how the hidden unit complexity of NQS depend on the local single-spin basis used. Exploiting the tensor network version of our RBM we construct an analytic NQS representation of the Affleck-Kennedy-Lieb-Tasaki (AKLT) state in the $xyz$ spin-1 basis using only $M=2N$ hidden units, compared to $M\sim O\left(N^2\right)$ required in the $S^z$ basis. Additional VMC calculations provide strong evidence that the AKLT state in fact possesses an exact compact NQS representation in the $xyz$ basis with only $M=N$ hidden units. These insights help to further unravel how to most effectively adapt the NQS framework for more complex quantum systems.

1. Introduction

The strongly-interacting quantum many-body problem is crucial to our understanding of many intriguing physical phenomena, but it is also inherently difficult to treat numerically owing to the exponential growth of the Hilbert space with system size. A commonly used approximate strategy is the variational method where a trial state, characterized by a tractable number of variational parameters, is optimized in energy. The effectiveness of this approach is highly dependent on the ansatz having an expressive form that can be systematically improved, to minimize bias, while also allowing relevant observables to be evaluated efficiently. Tensor networks have provided several examples of such ansatzes, with matrix product states (MPS) [1,2] displaying impressive accuracy in one-dimensional systems, along with projected entangled-pair states (PEPS) [1,3], tree tensor networks (TTN) [4,5], and multi-scale entanglement renormalization ansatz (MERA) [6] making two-dimensional systems accessible. Recently, artificial neural networks (ANNs) have emerged as another class of highly flexible variational ansatzes with many variants, such as restricted Boltzmann machines (RBM) [7], Deep Boltzmann Machines (DBM) [8,9,10], convolutional neural networks (CNN) [11,12,13,14,15], and feed-forward neural networks (FFNN) [16,17,18,19]. An important advantage of ANNs is that they are highly flexible and can be applied to any number of spatial dimensions, making them a powerful method for tackling the subtle physics seen in two-dimensional systems.

Although one of the simplest ANN variants, RBMs have seen widespread applications, including for open quantum systems [20,21,22,23], frustrated spin problems [24,25], quantum circuit simulation [26,27,28], and more. There are several reasons for their continued use. First, their simple structure allows for efficient sampling crucial for applying variational Monte Carlo (VMC) [29,30]. Second, RBMs are also a good candidate for a weakly biased ansatz, given that they are capable of exactly representing arbitrary states when their hidden unit number M scales exponentially with system size N. Third, RBMs are capable of representing states with volume-law entanglement [31], which further distinguishes them from tensor networks [32], despite their conceptual similarities [33,34,35]. Finally, there are also numerous classes of states with efficient exact RBM representations, including graph states [8]; spin Jastrow states, such as Laughlin states [33,36,37]; and general stabilizer states, such as the toric code [38,39,40,41], as well as more exotic hypergraph states and XS-stabilizer states [40]. Recently, we found that all but the last class listed, in fact, have RBM representations requiring $M<N$ hidden units [42], illustrating how even very modestly sized RBMs have significant representational power.

Despite their efficacy for spin-$\frac{1}{2}$ systems, the application of RBMs to systems with a local on-site dimension $d>2$, such as spin-1 or bosonic systems, has been limited with convolutional or feedforward neural networks generally being favored [16,17,43]. The typical approach in machine learning to handle models with multinomial or categorical variables is so-called “one-hot” or “unary” encoding [44]. Rather than representing a physical degree of freedom directly with one visible unit this approach encodes the possible local physical states into a set of binary visible units. While this approach leverages the power of binary or spin-$\frac{1}{2}$ RBMs, it multiplies the number of visible units by a factor d, significantly increasing the parameter count and complexity of the optimization. As a consequence, studies utilizing unary encoding so far, for example, on the Bose Hubbard model [45,46], have been limited to small system sizes. Thus, there is need to devise more efficient RBM constructions tailored for $d>2$ systems.

Progress has been made in this direction in recent work [47], where multivalued RBMs were applied directly to the one-dimensional spin-1 anti-ferromagnetic Heisenberg model (AFH) and substantially enhanced by incorporating a transformation to a coupled SU(2) symmetric basis. Complementary to this, here, we propose and study a direct generalization of the RBMs to spin-1 systems that retains key properties of spin-$\frac{1}{2}$ RBM with a minimal increase in variational parameters (as we will not be examining other network architectures, we will use RBM and NQS interchangeably in this paper). Specifically, the ability to describe arbitrary product states without hidden units, invariance of the parameterization to the values assigned to visible variables (labeling freedom), and equivalence to the tensor network formulation. This leads to the introduction of new quadratic bias and interaction weights in the RBM effective energy function. We demonstrate the effectiveness of the new formulation via VMC calculations for the spin-1 AFH model, where it is seen to deliver the same accuracy as unary encoding but with substantially fewer variational parameters. Additionally, we also investigate how the local single-spin basis affects the hidden unit complexity of a state by performing calculations in both the $S^z$ and $xyz$ spin-1 bases. For the AFH model with periodic boundary conditions, we find that the $S^z$ is more accurate. A useful advantage of our new spin-1 RBM formulation is that it permits tensor network based analytic constructions. Focusing on the paradigmatic Affleck-Kennedy-Lieb-Tasaki (AKLT) model, we give explicit exact NQS representations in both the $S^z$ and $xyz$ bases. Our $S^z$ basis NQS construction displays the expected [36] $M\sim O\left(N^2\right)$ scaling, while the simplification of the amplitude structure in the $xyz$ basis gives an NQS construction with $M=2N$ hidden units. By using VMC calculations, we find compelling evidence that the AKLT state, in fact, only requires $M=N$ hidden units to be represented exactly in the $xyz$ basis.

The structure of this paper is as follows. We briefly outline the VMC method applied to the many body problem in Section 2 and discuss desirable properties for variational ansatzes. Next, in Section 3, we discuss the RBM in its spin-$\frac{1}{2}$ form and analyze its key properties. In Section 4, we introduce a new generalization of NQS to spin-1 systems designed to mimic these properties and present VMC results for the AFH model. We then introduce analytic constructions for the AKLT model in Section 5, followed by VMC calculations in both the $S^z$ and $xyz$ basis. Finally, in Section 6, we conclude and discuss some open problems.

2. The Many-Body Problem and Variational Monte Carlo

2.1. Quantum Many-Body Problem

In this work, we will focus on physical systems composed of N localized spinful particles. Each particle is described by a vector of spin operators ${\widehat{\mathit{S}}}_j=({\widehat{S}}_j^x,{\widehat{S}}_j^y,{\widehat{S}}_j^z)$ for $j=1,2,\cdots ,N$. Typically, the eigenstates $\left|S_j\right.〉$ of ${\widehat{S}}_j^z$ are used as the local spin basis, from which the $S^z$ basis for the full system is constructed as $\left|\mathit{S}\right.〉=\left|S_1\right.〉\otimes \cdots \otimes \left|S_N\right.〉$, where $\mathit{S}=(S_1,S_2,\cdots ,S_N)$. Any many-body quantum state for these N spins can then be expanded in this basis as

$$\left|\mathsf{\Psi }\right.〉=\sum _{\mathit{S}}\mathsf{\Psi }\left(\mathit{S}\right)\left|\mathit{S}\right.〉,$$

(1)

via its complex amplitudes $\mathsf{\Psi }\left(\mathit{S}\right)$. In the spin-$\frac{1}{2}$ case, we have ${\widehat{S}}_j^{\alpha }=\frac{1}{2}{\widehat{\sigma }}_j^{\alpha }$ (taking $\hslash =1$ throughout) for $\alpha =\{x,y,z\}$ defined by the Pauli matrices, and $S_j\in \{+\frac{1}{2}\equiv \uparrow ,-\frac{1}{2}\equiv \downarrow \}$. In the spin-1 case, we have

$$\begin{array}{ccc}\hfill {\widehat{S}}_j^x& =& \frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 1\\ 0& 1& 0\end{array}\right],{\widehat{S}}_j^y=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0& -\mathrm{i}& 0\\ \mathrm{i}& 0& -\mathrm{i}\\ 0& \mathrm{i}& 0\end{array}\right],{\widehat{S}}_j^z=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right],\hfill \end{array}$$

with $S_j\in \{+1\equiv \Uparrow ,0,-1\equiv \Downarrow \}$.

A key challenge in many-body physics is to find the ground state $\left|{\mathsf{\Psi }}_0\right.〉$ of a system governed by an interacting Hamiltonian $\widehat{H}$. In the context of spin systems, $\widehat{H}$ will include terms that are products of the spin operators ${\widehat{\mathit{S}}}_j$ over two or more spins. Given that the expectation value of any observable $\widehat{A}$ for a general (unnormalized) state $\left|\mathsf{\Psi }\right.〉$ is

$${\langle A\rangle }_{\mathsf{\Psi }}=\frac{\left.〈\mathsf{\Psi }\right|\widehat{A}\left|\mathsf{\Psi }\right.〉}{〈\mathsf{\Psi }|\mathsf{\Psi }〉},$$

(2)

the variational approach reformulates the eigenvalue problem $\widehat{H}\left|{\mathsf{\Psi }}_0\right.〉=E_0\left|{\mathsf{\Psi }}_0\right.〉$ as the minimization of the energy $E={\langle \widehat{H}\rangle }_{\mathsf{\Psi }}$ over the exponentially many amplitudes $\mathsf{\Psi }\left(\mathit{S}\right)$. Performing this task exactly is only feasible for small systems $N\sim O\left(10\right)$ [48].

2.2. Variational Monte Carlo Method

A way to circumvent the “curse of dimensionality” is to instead to restrict the optimization over a specialized class of states $\left|{\mathsf{\Psi }}_{\mathit{p}}\right.〉$, dependent on parameters $\mathit{p}$ whose number $n_{\mathrm{params}}$ scales polynomially with N. The variational principle $E_0\le E_{{\mathit{p}}_0}={\min }_{\mathit{p}}{\langle \widehat{H}\rangle }_{{\mathsf{\Psi }}_{\mathit{p}}}$ then provides a route to finding the best approximation $\left|{\mathsf{\Psi }}_{{\mathit{p}}_0}\right.〉$ within the ansatz for $\left|{\mathsf{\Psi }}_0\right.〉$. The flexibility and utility of variational ansatz are greatly enhanced if instead of computing expectation values ${\langle A\rangle }_{{\mathsf{\Psi }}_{\mathit{p}}}$ exactly we evaluate them approximately using Monte Carlo sampling. This approach is called variational Monte Carlo (VMC) and is described in detail in Appendix A. A key feature of VMC is that only ratios of amplitudes for an ansatz ${\mathsf{\Psi }}_{\mathit{p}}\left(\mathit{S}\right)/{\mathsf{\Psi }}_{\mathit{p}}\left({\mathit{S}}^{\prime }\right)$ between different spin states $\left|\mathit{S}\right.〉$ and $\left|{\mathit{S}}^{\prime }\right.〉$ are needed within the algorithm allowing us to ignore the normalization of quantum states in this work. To be efficient, thus, we require that these amplitude ratios for our ansatz can be evaluated with an effort scaling polynomially with N.

All numerical calculations presented in this paper employ the powerful stochastic reconfiguration method [30,49] for optimizing the parameters $\mathit{p}$. While ansatzes with a larger number of parameters may describe more varied amplitude structure, in principle allowing greater accuracy, the effort to optimize them can scale as $O\left(n_{\mathrm{params}}^3\right)$. It is, therefore, crucial that the functional form of any candidate ansatz $\left|{\mathsf{\Psi }}_{\mathit{p}}\right.〉$ is judicious in its use of variational parameters to avoid excessive redundancies.

3. Neural-Network Quantum States in Spin-1/2

In this section, we review NQS in terms of RBMs, discussing their form from both an energy function and tensor network perspective. In doing so, we highlight key properties of NQS that are desirable for expressiveness and the ability to represent important families of states.

3.1. Restricted Boltzmann Machine Approach

Restricted Boltzmann machines consist of two sets of units, Nvisible units representing the physical system, and Mhidden units to be marginalized out. The units are characterized by a Boltzmann-like combined probability distribution $p(\mathit{v},\mathit{h})=exp(-{\mathcal{E}}_{\mathbf{\lambda }})$ with an effective “energy” function

$${\mathcal{E}}_{\mathbf{\lambda }}(\mathit{v},\mathit{h})=-\sum _{j=1}^N{a}_j{v}_j-\sum _{i=1}^M{b}_i{h}_i-\sum _{i=1}^M\sum _{j=1}^N{w}_{ij}{h}_i{v}_j,$$

(3)

where $\mathbf{\lambda }=\{\mathit{a},\mathit{b},\mathit{w}\}$ is the set of $MN+M+N$ complex parameters consisting of N visible biases $\mathit{a}=(a_1,\cdots ,a_N)$, M hidden biases $\mathit{b}=(b_1,\cdots ,b_M)$, and $M\times N$ hidden-visible interaction weights $\mathit{w}=[w_{11},\cdots ,w_{1N},\cdots ,w_{MN}]$. While $p(\mathit{v},\mathit{h})$ is Boltzmann once the hidden units are traced out, we are left with a more complex marginal distribution for the visible units, from which the NQS amplitudes are derived [7]:

$${\mathsf{\Psi }}_{\mathbf{\lambda }}(\mathit{S}\mapsto \mathit{v})=\sum _{\mathit{h}}exp\left[-{\mathcal{E}}_{\mathbf{\lambda }}(\mathit{v},\mathit{h})\right].$$

(4)

The RBM architecture is shown diagrammatically in Figure 1.

Typically, visible $v_j$ and hidden $h_i$ unit variables are taken as two-valued. The pair of unique values $\{\mu ,\nu \}$ taken by units within the energy function Equation (3) can be freely chosen, and they need not coincide with the eigenvalues $S_j\in \{+\frac{1}{2},-\frac{1}{2}\}$ of the local operator ${\widehat{S}}^z$ chosen to define the physical basis. In Equation (4), we emphasize this labeling freedom by explicitly introducing a mapping between physical and visible configuration $\mathit{S}\mapsto \mathit{v}$. Commonly, with NQS, an implicit choice $\mathit{v}=\mathit{S}$ is made, nullifying the need for this distinction, but it will prove to be useful here. Canonical choices $v_j$ are number-like $\{0,1\}$ or Ising-like $\{+1,-1\}$ values. We will adopt Ising-like visible variables, in which case it follows that an arbitrary visible unit taking values

$$S_j=\left\{\begin{array}{ccc}\uparrow & \mapsto & \hfill \mu \\ \downarrow & \mapsto & \hfill \nu \end{array}\right\}={\tilde{v}}_j$$

(5)

is generated by a shift and rescaling of $v_j\in \{+1,-1\}$ as ${\tilde{v}}_j=\frac{1}{2}(\mu +\nu )+\frac{1}{2}(\mu -\nu )v_j$. Since this is a linear transformation, it is entirely accommodated by redefining the energy function ${\mathcal{E}}_{\mathbf{\lambda }}(\mathit{v},\mathit{h})$ parameters $\mathbf{\lambda }$ without changing its functional form. We shall see shortly that labeling freedom is a crucial ingredient for generalizing RBMs to spin-1 system using higher-dimensional visible units in Section 4. Tracing out Ising-like hidden variables results in the amplitude expansion [7]

$$\begin{array}{ccc}\hfill {\mathsf{\Psi }}_{\mathbf{\lambda }}(\mathit{S}\mapsto \mathit{v})& =& \prod _{j=1}^N{e}^{a_j{v}_j}\prod _{i=1}^{M}2cosh\left(b_i+\sum _{j=1}^N{w}_{ij}{v}_j\right),\hfill \end{array}$$

(6)

that is commonly used and numerically convenient.

A key property of spin-$\frac{1}{2}$ RBMs is that they can represent any generic product state of the form

$$\left(c_{+}^{\left(1\right)}\left|\uparrow \right.〉+c_{-}^{\left(1\right)}\left|\downarrow \right.〉\right)\otimes \cdots \otimes \left(c_{+}^{\left(N\right)}\left|\uparrow \right.〉+c_{-}^{\left(N\right)}\left|\downarrow \right.〉\right),$$

(7)

without hidden units by simply setting their visible biases to

$$a_j=\frac{ln\left(c_{+}^{\left(j\right)}\right)-ln\left(c_{-}^{\left(j\right)}\right)}{\mu -\nu }.$$

(8)

Hidden units are, thus, necessary to describe entangled states; however, unlike tensor networks, there is no direct quantitative relation between them. Indeed, Ref. [31] reports that, for random NQS, adding hidden units actually decreases the amount of entanglement in a state. Nevertheless, increasing the number M of hidden units increases the expressiveness of the NQS, allowing more complex correlations within ${\mathsf{\Psi }}_{\mathbf{\lambda }}(\mathit{S}\mapsto \mathit{v})$ to be encoded. Formally, the NQS ansatz can represent any arbitrary state exactly in the limit of $M\to 2^N$[50]. More usefully, there are important classes of states that possess highly accurate approximate or even exact NQS representations with an efficient scaling $M\sim \mathrm{poly}\left(N\right)$  [33,36,38,39,40,41]. A particularly tractable NQS hidden unit complexity is typified by the following: efinition

Definition 1

(Compact NQS). States with an exact NQS representation where $M\le N$ will be denoted as compact.

As we demonstrated in ref. [42], important classes of state, including Jastrow, graph, and stabilizer states, all have compact NQS representations.

3.2. Tensor Network Approach

An alternative formulation of NQS views them as a tensor network. From this perspective, the amplitudes $\mathsf{\Psi }\left(\mathit{S}\right)$ are recast as elements of an order-N tensor ${\mathsf{\Psi }}_{S_1{S}_2\cdots S_N}$, and this structureless tensor can be decomposed into a set of lower order tensors contracted together [1,2]. Here, we will make repeated use of tensor network diagrams that form an important analytical tool. Here, generic tensors of any order are represented as shapes, most often a circle ∘ color shaded to guide the eye, with protruding legs for each index they possess. Contraction of tensors is then represented by the joining of legs via graphical equations like this:

(9)

depicting $C_{abxz}={\sum }_{\alpha }{A}_{ab\alpha }{B}_{x\alpha z}$. When there is a symbol inside a shape, it represents tensors with fixed elements or some specific structure. For example, the order-1 tensors with ↑ or ↓ symbols are the representation of the spin-$\frac{1}{2}$$S^z$ basis states $\left|\uparrow \right.〉=\left(\genfrac{}{}{0pt}{}{1}{0}\right)$ and $\left|\downarrow \right.〉=\left(\genfrac{}{}{0pt}{}{0}{1}\right)$, while the tensor shown in Figure 2a represents $\left|+\right.〉=\left|\uparrow \right.〉+\left|\downarrow \right.〉$. A special tensor we will use frequently is the COPY tensor [33,51], denoted by a dot •, which is the multidimensional generalization of the $2\times 2$ identity matrix. Its order-3 variant is shown in Figure 2b. It possesses a number of important properties. First, contracting any COPY tensor leg with $\left|+\right.〉$ reduces the order of the COPY tensor by 1, as shown in Figure 2c. Second, contracting any COPY tensor leg with an $S^z$ basis state factorizes it by duplicating the basis vector across all legs, as shown in Figure 2d, motivating its name. This property makes COPY tensors a useful glue for constructing sampleable tensor networks. An order-N COPY tensor can itself be viewed as a representation of a GHZ state, shown in Figure 2e. Finally, Figure 2f shows a generic diagonal matrix attached to a leg of a COPY tensor can be commuted across to any of the COPY dot’s other legs.

With these concepts in place, the bipartite graph structure of RBM graph readily translates to a tensor network:

(10)

in which the vertices representing the visible and hidden units are replaced by COPY tensors, and order-2 tensors, or $2\times 2$ coupling matrices ${\mathbf{C}}^{\left(ij\right)}$, are introduced on each edge of the graph. Based on Figure 2e, thus, we can view each hidden unit in an NQS as contributing within an amplitude-wise product a locally deformed GHZ state. The NQS tensor network readily displays the key properties of the spin-$\frac{1}{2}$ RBM. They can trivially represent product states in Equation (7) by using one hidden unit with rank-1 coupling matrices $\mathit{P}$ containing the coefficients $\{c_{+}^{\left(j\right)},c_{-}^{\left(j\right)}\}$ that subsequently factorizes out as

(11)

Moreover, since the elements of coupling matrices ${\mathbf{C}}^{\left(ij\right)}$ have no functional dependence on any visible variables, they manifestly display labeling freedom. Crucially, while the NQS tensor network appears to have $4MN$ variational parameters, it exhibits significant gauge freedom due to the ability to reshuffle diagonal matrices across COPY tensors, as shown in Figure 2g. Consequently, most of the elements in the coupling matrices can be extracted and combined at the COPY tensors, reducing the number of independent parameters to $N+M+NM$, identical to the RBM formulation [42]. This equivalence allows either the NQS tensor network or RBM approach to be used interchangeably, and any direct generalization of RBMs to spin-1 should retain this feature.

4. Generalization to Spin-1 Systems

A tentative definition of a spin-1 RBM can be made by simply using the three-valued visible variables $v_j$ with a natural Ising-like physical to visible mapping,

$$S_j=\left\{\begin{array}{ccc}\Uparrow & \mapsto & \hfill +1\\ 0& \mapsto & \hfill 0\\ \Downarrow & \mapsto & \hfill -1\end{array}\right\}=v_j,$$

(12)

in the standard linear energy function Equation (3). By retaining two-valued Ising-like hidden units $h_i\in \{+1,-1\}$, the amplitudes ${\mathsf{\Psi }}_{\mathbf{\lambda }}(\mathit{S}\mapsto \mathit{v})={\sum }_{\mathit{h}}exp\left[{\mathcal{E}}_{\mathbf{\lambda }}(\mathit{v},\mathit{h})\right]$ continue to be given by Equation (6), but now admitting three-valued visible variables.

This spin-1 generalization of the RBM lacks the key properties of the spin-$\frac{1}{2}$ variant. First, ${\mathsf{\Psi }}_{\mathbf{\lambda }}(\mathit{S}\mapsto \mathit{v})$ cannot easily describe a generic product spin state

$$\left(c_{+}^{\left(1\right)}\left|\Uparrow \right.〉+c_0^{\left(1\right)}\left|0\right.〉+c_{-}^{\left(1\right)}\left|\Downarrow \right.〉\right)\otimes \cdots \otimes \left(c_{+}^{\left(N\right)}\left|\Uparrow \right.〉+c_0^{\left(N\right)}\left|0\right.〉+c_{-}^{\left(N\right)}\left|\Downarrow \right.〉\right),$$

(13)

with coefficients $c_{\alpha }^{\left(j\right)}$, since the visible bias term $exp\left(a_j{v}_j\right)$ cannot discriminate each of the local spin-1 states. This points to issues of expressiveness since we are using the same set of $MN+M+N$ parameters $\{\mathit{a},\mathit{b},\mathit{w}\}$ to describe a bigger state space. Second, the energy function Equation (3) does not possess labeling freedom for the visible units values. An arbitrary physical to visible mapping,

$$S_j=\left\{\begin{array}{ccc}\Uparrow & \mapsto & \hfill \mu \\ 0& \mapsto & \hfill \nu \\ \Downarrow & \mapsto & \hfill \omega \end{array}\right\}={\tilde{v}}_j,$$

(14)

is generated from an Ising-like variable $v_j=\{+1,0,-1\}$ via the quadratic transformation ${\tilde{v}}_j=\nu +\frac{1}{2}(\mu -\omega )v_j+\frac{1}{2}(\mu +\omega -2\nu )v_j^2$. Consequently, transforming visible variables in this way cannot be accommodated in ${\mathcal{E}}_{\mathbf{\lambda }}(\mathit{v},\mathit{h})$ by changing only the parameters $\{\mathit{a},\mathit{b},\mathit{w}\}$. One way to avoid these issues is to use “one-hot” or “unary” encoding which has been successfully applied to both spin-1 [40] and bosonic systems [45,46].

4.1. Unary Encoding Approach

Unary encoding applies the spin-$\frac{1}{2}$ RBM formalism of Equation (4) to a larger state space of each spin-1 by mapping the physical system into one comprising a larger number of spin-$\frac{1}{2}$ particles. Specifically, the local $S^z$ basis of each spin-1 is mapped on to three spin-$\frac{1}{2}$’s as

$$\left|\Uparrow \right.〉=\left|\uparrow \uparrow \downarrow \right.〉,\left|0\right.〉=\left|\uparrow \downarrow \uparrow \right.〉,\left|\Downarrow \right.〉=\left|\downarrow \uparrow \uparrow \right.〉.$$

(15)

Physical states of the spin-1 system are now contained in the unary encoded subspace where only a single-excitation occurs within any three-site unit cell which facilitates the efficient projection of the representation. This naturally generalizes to a d-dimensional local state space.

By using this mapping, a spin-$\frac{1}{2}$ RBM can be applied inheriting its labeling freedom. Owing to the single-excitation projection, product states in Equation (13) are readily described by the visible biases $a_{j\mathrm{a}},a_{j\mathrm{b}}$ and $a_{j\mathrm{c}}$ associated with the three spin-$\frac{1}{2}$’s {a,b,c} encoding a given spin-1. In general, for a d-dimensional local state space, unary encoding accounts for the enlarged state space by increasing the number of variational parameters to $n_{\mathrm{params}}=M+dN+dNM$ compared to the naive spin-1 RBM. However, there are some obvious deficiencies of this approach. First, unary encoding appears to unnecessarily enlarge the parameter count, as evidenced by the fact that it increases it even for $d=2$. This will significantly increase the computational cost. Second, splitting the interaction weights w across a unary cell makes the interpretation of any individual hidden unit’s contribution to the physical state rather opaque.

4.2. Defining a Spin-1 RBM and Tensor Network

These deficiencies show that a more general RBM energy function is needed that is sensitive to the three-valued nature of the visible units through the inclusion of terms involving the square of visible variables (a similar approach is likely to have been used in Ref. [47] already, although it was not explicitly stated). This motivates the following definition: efinition

Definition 2

(spin-1 NQS). We introduce a direct spin-1 RBM as the ansatz with amplitudes ${\mathsf{\Psi }}_{\mathsf{\Lambda }}(\mathit{S}\mapsto \mathit{v})={\sum }_{\mathit{h}}exp\left[{\mathcal{E}}_{\mathsf{\Lambda }}(\mathit{v},\mathit{h})\right]$ via the energy function

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathsf{\Lambda }}(\mathit{v},\mathit{h})=\sum _{j=1}^N{a}_j{v}_j+\sum _{j=1}^N{A}_j{v}_j^2+\sum _{i=1}^M{b}_i{h}_i+\sum _{i=1}^M\sum _{j=1}^N{w}_{ij}{h}_i{v}_j+\sum _{i=1}^M\sum _{j=1}^N{W}_{ij}{h}_i{v}_j^2,\end{array}$$

(16)

defined by the parameters $\mathsf{\Lambda }=\{\mathit{a},\mathit{b},\mathit{w},\mathit{W},\mathit{A}\}$.

The new contributions to a spin-1 RBM are $\mathit{W}$, an $M\times N$-dimensional matrix of quadratic interactions, and $\mathit{A}$, an N-dimensional vector of quadratic visible biases. There are now $2MN+2N+M$ complex parameters in total. Tracing out the two-valued Ising-like hidden units gives amplitudes

$$\begin{array}{ccc}\hfill {\mathsf{\Psi }}_{\mathsf{\Lambda }}(\mathit{S}\mapsto \mathit{v})& =& \prod _{i=1}^N{e}^{a_i{v}_i+A_i{v}_i^2}\prod _{j=1}^{M}2cosh\left(b_j+\sum _{i=1}^N{w}_{ij}{v}_i+\sum _{i=1}^N{W}_{ij}{v}_i^2\right).\hfill \end{array}$$

(17)

The inclusion of a quadratic visible bias now allows any product state in Equation (7) to be described without hidden units by setting

$$a_j=log\left(\sqrt{c_{+}^{\left(j\right)}{c}_{-}^{\left(j\right)}}/c_0^{\left(j\right)}\right)\mathrm{and}{A}_j=log\left(\sqrt{c_{+}^{\left(j\right)}/c_{-}^{\left(j\right)}}\right),$$

while the quadratic interaction term ensures labeling freedom for the visible variable.

A strong justification for Equation (16) being the appropriate spin-1 generalization of RBMs is its relation to an NQS tensor network for spin-1. To handle spin-1 systems, we introduce a COPY tensor with three-dimensional legs copying the $S^z$ basis states $\{\left|\Uparrow \right.〉,\left|0\right.〉,\left|\Downarrow \right.〉\}$. Its properties are summarized in Figure 3a–f and are straightforward generalizations of the spin-$\frac{1}{2}$ case given in Figure 2a–f and discussed in Section 3.2. The major difference is that we now distinguish three-dimensional legs with ‘=’ lines instead of ‘−’. As before, the NQS tensor network follows from the conversion of the RBM graph in Figure 1, except that visible vertices are now replaced by three-dimensional COPY tensors giving:

(18)

Owing to the mixed dimensionality of the COPY tensors in this network, it now requires a rectangular $2\times 3$ coupling matrix ${\mathbf{C}}^{\left(ij\right)}$ between the i-th hidden and j-th visible unit. For the i-th hidden unit, its set of coupling matrices ${{\rm Y}}^{\left(i\right)}=\{{\mathbf{C}}^{\left(i1\right)},{\mathbf{C}}^{\left(i2\right)},\cdots ,{\mathbf{C}}^{\left(iN\right)}\}$ can be explicitly tabulated as

$$\begin{array}{ccccc}\hfill & \stackrel{v_1}{\overbrace{+10-1}}& \stackrel{v_2}{\overbrace{+10-1}}& \cdots & \stackrel{v_N}{\overbrace{+10-1}}\\ \\ \hfill h_i\left\{\begin{array}{c}+1\\ -1\end{array}\right.& \left[\begin{array}{ccc}{C}_{++}^{\left(i1\right)}& C_{+0}^{\left(i1\right)}& C_{+-}^{\left(i1\right)}\\ C_{-+}^{\left(i1\right)}& C_{-0}^{\left(i1\right)}& C_{--}^{\left(i1\right)}\end{array}\right]& \left[\begin{array}{ccc}{C}_{++}^{\left(i2\right)}& C_{+0}^{\left(i2\right)}& C_{+-}^{\left(i2\right)}\\ C_{-+}^{\left(i2\right)}& C_{-0}^{\left(i2\right)}& C_{--}^{\left(i2\right)}\end{array}\right]& \cdots & \left[\begin{array}{ccc}{C}_{++}^{\left(iN\right)}& C_{+0}^{\left(iN\right)}& C_{+-}^{\left(iN\right)}\\ C_{-+}^{\left(iN\right)}& C_{-0}^{\left(iN\right)}& C_{--}^{\left(iN\right)}\end{array}\right]\end{array}.$$

The amplitudes ${{\rm Y}}^{\left(i\right)}\left(\mathit{v}\right)$ of this hidden unit’s correlator then follow by summing the product of coupling matrix elements selected by $\mathit{v}$ along each row, giving

$${{\rm Y}}^{\left(i\right)}\left(\mathit{v}\right)=\sum _{h_i=\pm 1}^1\prod _{j=1}^N{C}_{h_j{v}_i}^{\left(ij\right)}=C_{+v_1}^{\left(i1\right)}{C}_{+v_2}^{\left(i2\right)}\cdots C_{+v_N}^{\left(iN\right)}+C_{-v_1}^{\left(i1\right)}{C}_{-v_2}^{\left(i2\right)}\cdots C_{-v_N}^{\left(iN\right)}.$$

(19)

The amplitudes of NQS tensor network are then the product of each of these hidden unit correlators

$$\begin{array}{ccc}\hfill {\mathsf{\Psi }}_{\mathrm{NQS}}(\mathit{S}\to \mathit{v})& =& \prod _{i=1}^M{{\rm Y}}^{\left(i\right)}\left(\mathit{v}\right),\hfill \end{array}$$

(20)

which, like an RBM, can be exactly and efficiently sampled.

The spin-1 NQS tensor network appears to have $5MN$ complex variational parameters; however, again, gauge freedom allows the shuffling of diagonal matrices (of an appropriate dimension) through a COPY tensors reducing this. Specifically, its equivalence to the generalized spin-1 RBM proposed in Equation (16) is made using the coupling matrix decomposition

$$\begin{array}{ccc}\hfill \left[\begin{array}{ccc}{C}_{++}^{\left(ij\right)}& C_{+0}^{\left(ij\right)}& C_{+-}^{\left(ij\right)}\\ C_{-+}^{\left(ij\right)}& C_{-0}^{\left(ij\right)}& C_{--}^{\left(ij\right)}\end{array}\right]& =& e^c\left[\begin{array}{cc}{e}^{{\tilde{b}}_{ij}}& 0\\ 0& e^{-{\tilde{b}}_{ij}}\end{array}\right]\left[\begin{array}{ccc}{e}^{w_{ij}+W_{ij}}& 1& e^{-w_{ij}+W_{ij}}\\ e^{-w_{ij}-W_{ij}}& 1& e^{w_{ij}-W_{ij}}\end{array}\right]\hfill \\ & & \times \left[\begin{array}{ccc}{e}^{{\tilde{a}}_{ij}+{\tilde{A}}_{ij}}& 0& 0\\ 0& 1& 0\\ 0& 0& e^{-{\tilde{a}}_{ij}+{\tilde{A}}_{ij}}\end{array}\right].\hfill \end{array}$$

The solution for the weights $w_{ij},W_{ij}$ and partial biases ${\tilde{a}}_{ij},{\tilde{A}}_{ij},{\tilde{b}}_{ij}$ is outlined in Appendix B, from which the full RBM biases are formed as $a_j={\sum }_{i=1}^M{\tilde{a}}_{ij}$, $A_j={\sum }_{i=1}^M{\tilde{A}}_{ij}$ and $b_i={\sum }_{j=1}^N{\tilde{b}}_{ij}$. The spin-1 NQS tensor network, thus, reduces to $n_{\mathrm{params}}=2MN+2N+M$ complex parameters.

This correspondence between the spin-1 RBM and the spin-1 NQS tensor network highlights an advantage over unary encoding. Specifically, coupling matrices provide an intuitive tool for engineering the correlations and structures that a given hidden unit imprints on the amplitudes of an NQS. A trivial case is when all the elements of the j-th coupling matrix are 1’s, denoted generically as $\mathbf{I}$, equivalent to the hidden unit being disconnected from that visible unit. A more complex example with conditional correlations is a hidden unit with coupling matrices

$$\{\mathbf{I},\cdots ,\mathbf{I},\stackrel{k}{\overbrace{{\mathbf{C}}_{\mathrm{c}\Uparrow }}},{\mathbf{C}}_{-0\Downarrow },\cdots ,{\mathbf{C}}_{-0\Downarrow }\},$$

built from

$${\mathbf{C}}_{\mathrm{c}\Uparrow }=\left[\begin{array}{ccc}0& 1& 1\\ 1& 0& 0\end{array}\right]\mathrm{and}{\mathbf{C}}_{-0\Downarrow }=\left[\begin{array}{ccc}1& \hfill 1& \hfill 1\\ 1& \hfill -1& \hfill -1\end{array}\right].$$

This hidden unit generates a correlator ${\rm Y}\left(\mathit{v}\right)={\delta }_{v_k,\Downarrow }+{\delta }_{v_k,\Uparrow }{(-1)}^{\mathcal{C}\left(\mathit{v}\right)}$, where a factor ${(-1)}^{\mathcal{C}\left(\mathit{v}\right)}$ is introduced conditional on the kth spin being in the state $\left|\Uparrow \right.〉$, with $\mathcal{C}\left(\mathit{v}\right)$ being the number of $\left|0\right.〉$ and $\left|\Downarrow \right.〉$ states in the configuration $\mathit{v}$ between spins $k+1$ and N. Similar types of hidden units will be used extensively in Section 5.1 to construct an exact representation of a state.

4.3. Projection of Unary Encoding into a Spin-1 RBM

The tensor network formalism provides further evidence that unary encoding from Section 4.1 is an over-parameterization of RBMs for spin-1 systems with $\delta n_{\mathrm{params}}=N(1+M)$ redundant parameters. Unary projection is implemented by an order-4 tensor $\mathcal{U}$ obeying:

(21)

By directly applying this projector to the spin-$\frac{1}{2}$ NQS tensor network for a unary encoded state and performing a graphical rewrite, we obtain the spin-1 NQS introduced in Equation (18). Figure 4 summarizes the crucial manipulations required.

In Figure 4a, some representative examples of contractions between the projection tensor $\mathcal{U}$ and hidden units are shown. There are three steps to rewriting the network. The first step, shown in Figure 4b, essentially pulls $\mathcal{U}$ through the three unary two-dimensional COPY tensors, leaving behind a single three-dimensional COPY tensor representing the physical spin-1. Two important cases are shown in the example in Figure 4b. A hidden unit may have connections to each of the unary spin-$\frac{1}{2}$’s, where upon they get bundled up by the $\mathcal{U}$ tensor. A hidden unit may connect to only a subset of the unary spin-$\frac{1}{2}$’s, which is handled by plugging the unused legs of $\mathcal{U}$ tensor with $\left|+\right.〉$. If a hidden unit couples to more than one of the unary spins, then the second step, shown in Figure 4c, involves splitting the hidden unit’s two-dimensional COPY tensor to separate those connections. The final step is then to contract the split COPY tensor, coupling matrices and the projection $\mathcal{U}$ to form a rectangular $2\times 3$ coupling matrix, as depicted in Figure 4d. If a hidden unit has connections exclusively within the unary spin-$\frac{1}{2}$’s, then it becomes an entirely local visible bias contribution in the spin-1 NQS.

4.4. Change of Local Spin Basis

For tensor network representations, such as MPS or PEPS, the complexity (internal bond dimension) of a given state’s description is rooted in its entanglement structure. As such changing the local basis of the spins used in a calculation has no effect on this complexity. Moreover, transforming a representation from one basis to another is accomplished by simply admixing the local tensors. Within VMC, a change of local spin basis leaves the locality and the sparsity of the Hamiltonian essentially unchanged. However, it is pivotal to the method that the amplitudes ${\mathsf{\Psi }}_{\mathit{p}}\left(\mathit{S}\right)$ of whatever ansatz is used can be efficiently evaluated in this new basis. This is not generally true of NQS since their sampleability is intimately tied to the basis that factorizes the COPY tensors they are built from.

To understand how NQS behave, consider a representation of some state $\left|\mathsf{\Psi }\right.〉$, for instance, of four spin-1’s sampleable in the $S^z$ basis

(22)

Now, suppose we transform the local $S^z$ basis $\left|S\right.〉$ to a new basis $\chi$ via a unitary $\left|{\chi }_S\right.〉={\widehat{\mathcal{B}}}^{†}\left|S\right.〉$. Formally, from a tensor network perspective, we find the $\chi$ basis NQS representation by sandwiching $𝟙=\widehat{\mathcal{B}}{\widehat{\mathcal{B}}}^{†}$ on each physical leg and computing the $S^z$ basis NQS representation of ${\left({\widehat{\mathcal{B}}}^{†}\right)}^{\otimes N}\left|\mathsf{\Psi }\right.〉$. However, currently, there is no known procedure for updating exactly an NQS after the application of an arbitrary single-spin unitary, even allowing for an increased number of hidden units. While we are guaranteed that an NQS representation of ${\widehat{\mathcal{B}}}^{\otimes N}\left|\mathsf{\Psi }\right.〉$ exists, as illustrated here schematically

(23)

there is no guarantee it will be efficient, even if the representation of $\left|\mathsf{\Psi }\right.〉$ was originally. An example of such a catastrophic loss of NQS efficiency has been presented by Gao and Duan [8]. They show that, so long as the polynomial hierarchy in computational complexity theory does not collapse, a two-dimensional cluster state, which has an efficient NQS representation, has no efficient NQS representation after a specific layer of translation-invariant single-spin unitaries are applied. Thus, NQS complexity depends non-trivially on the local spin basis used.

Motivated by this, we will consider NQS calculations in two different local spin-1 bases to examine how the complexity varies, specifically, the standard $S^z$ basis and the $xyz$ basis defined as

$$\begin{array}{c}\hfill \left|x\right.〉=\frac{1}{\sqrt{2}}\left(\left|\Uparrow \right.〉-\left|\Downarrow \right.〉\right),\left|y\right.〉=\frac{\mathrm{i}}{\sqrt{2}}\left(\left|\Uparrow \right.〉+\left|\Downarrow \right.〉\right),\left|z\right.〉=-\left|0\right.〉.\end{array}$$

In the $xyz$ basis, the individual spin-1 operators all acquire the same off-diagonal form

$$\begin{array}{ccc}\hfill {\widehat{S}}_j^x& =& \left[\begin{array}{ccc}0& 0& 0\\ 0& 0& \mathrm{i}\\ 0& -\mathrm{i}& 0\end{array}\right],{\widehat{S}}_j^y=\left[\begin{array}{ccc}0& 0& \mathrm{i}\\ 0& 0& 0\\ -\mathrm{i}& 0& 0\end{array}\right],{\widehat{S}}_j^z=\left[\begin{array}{ccc}0& \mathrm{i}& 0\\ -\mathrm{i}& 0& 0\\ 0& 0& 0\end{array}\right],\hfill \end{array}$$

as a consequence of them all contributing one eigenstate to the basis. On a practical level, using the $xyz$ basis for the spin-1 RBM in Equation (17) simply requires replacing $\mathit{S}$ with $\mathbf{\alpha }$ and a physical-visible mapping from $\mathbf{\alpha }\mapsto \mathit{v}$, such as

$${\alpha }_j=\left\{\begin{array}{ccc}x& \mapsto & \hfill +1\\ y& \mapsto & \hfill 0\\ z& \mapsto & \hfill -1\end{array}\right\}=v_j.$$

(24)

Equivalently, the visible unit COPY tensors for this NQS tensor network can be considered to be rotated into this $xyz$ basis.

4.5. Numerical Example—Spin-1 Anti-Ferromagnetic Heisenberg Model

To confirm the effectiveness of our spin-1 NQS, we performed VMC calculations to reach the ground state of the well-known anti-ferromagnetic Heisenberg (AFH) model in one dimension. The Hamiltonian is given by

$$\widehat{\mathcal{H}}=J\sum _i{\widehat{\mathit{S}}}_i·{\widehat{\mathit{S}}}_{i+1},$$

(25)

where $J>0$ is the magnetic interaction strength and with periodic boundary conditions $N+1\equiv 1$. We focus on small systems allowing direct comparison to the ground state calculated from exact diagonalization via the overlap $\mathcal{O}=|〈{\mathsf{\Psi }}_{\mathrm{ED}}|{\mathsf{\Psi }}_{\mathsf{\Lambda }}〉{|}^2$. Additionally, we performed NQS optimizations in both the $S^z$ and $xyz$ bases to compare any differences in performance. This basis change alters the fixed quantum numbers of the system. In particular, in the $S^z$ basis, the AFH model preserves the total $S^z$ projection of the system, and its ground state lies in the ${\sum }_{j=1}^N{S}_j^z=0$ sector of the full Hilbert space, while, in the $xyz$ basis, the AFH model preserves the “parity” of the total $x,y,z$ spin populations in the system, and, for even N, the ground state lies in the subspace of configurations, where there is an even number of each basis state. As is common for VMC calculations, we only select configurations from the relevant subspace during sampling. We perform the calculations with hidden unit numbers $M=[1,2,\cdots ,2N]$, initializing $M=1$ with random small complex parameters. For each successive calculation, we need the NQS with the parameters for $M-1$ and initialize the Mth hidden unit with random small parameters, gradually increasing the size of the network in a sequential manner. To check the robustness against initialization bias of the qualitative features we have discussed, we rerun optimization sequences 5 to 10 times and present the best results here.

In Figure 5, we show how the accuracy of spin-1 NQS and unary encoding representations improve with an increasing number of hidden units M for both the $S^z$ and $xyz$ bases plotted in terms of the variational parameter count. The spin-1 NQS achieves a superior accuracy to unary encoding for a similar $n_{\mathrm{params}}$. The inset of Figure 5a shows the collapse of the same $S^z$ data plotted against M, indicating that the spin-1 NQS and unary encoding have in fact located the same solution for a given M. However, by using $\delta n_{\mathrm{params}}=N(1+M)$ less parameters, the spin-1 NQS is considerably more efficient to optimize, especially noting that $\delta n_{\mathrm{params}}$ scales with both system size and hidden unit number (for example, consider that the $N=12,M=12$ spin-1 NQS has a comparable parameter count to an $N=12,M=8$ unary NQS). In Figure 5b, we observe a noticeable drop in accuracy for both NQS variants in the $xyz$ basis compared to the $S^z$ basis. This suggests that the AFH ground state amplitude structure with periodic boundaries is inherently more complicated in the $xyz$ basis regardless of encoding. Moreover, this confirms that hidden unit number M of an NQS is basis-dependent quantity and cannot be used as a proxy of the entanglement.

5. Revisiting the AKLT Model

We now move on to benchmark our spin-1 NQS against the analytically solvable AKLT model [52], which is a spin-1 chain governed by a bilinear-biquadratic SU(2)-isotropic Heisenberg Hamiltonian of the form

$${\widehat{H}}_{\mathrm{AKLT}}=\sum _{j=1}^N\left[{\widehat{\mathit{S}}}_j·{\widehat{\mathit{S}}}_{j+1}+\beta {\left({\widehat{\mathit{S}}}_j·{\widehat{\mathit{S}}}_{j+1}\right)}^2\right]+\frac{2}{3},$$

with periodic boundary conditions $N+1\equiv 1$. It has special significance since it was the first solvable spin-1 chain model that exhibits the ‘Haldane gap’ [53]. The AKLT state $\left|{\mathsf{\Psi }}_{\mathrm{AKLT}}\right.〉$ is the ground state of ${\widehat{H}}_{\mathrm{AKLT}}$ at the AKLT point $\beta =\frac{1}{3}$, and has an energy of exactly zero.

As is well-known, $\left|{\mathsf{\Psi }}_{\mathrm{AKLT}}\right.〉$ has a special structure of correlations which are related to a valence bond solid. Specifically, each spin-1 is envisaged as being a pair of spin-$\frac{1}{2}$ particles that are correspondingly entangled in a singlet state with a partner spin-$\frac{1}{2}$ in the nearest neighboring spin-1 on the chain. The AKLT state is then the projection $\mathbb{P}$ of the local pairs of spin-$\frac{1}{2}$ particles into the triplet subspace, as depicted in Figure 6. This also leads to $\left|{\mathsf{\Psi }}_{\mathrm{AKLT}}\right.〉$ possessing a very compact MPS representation with matrices

$$\begin{array}{ccc}\hfill {\mathbf{A}}^{+}& =& \frac{1}{\sqrt{2}}{\widehat{\sigma }}^{+},{\mathbf{A}}^0=-\frac{1}{2}{\widehat{\sigma }}^z,{\mathbf{A}}^{-}=-\frac{1}{\sqrt{2}}{\widehat{\sigma }}^{-},\hfill \end{array}$$

where ${\widehat{\sigma }}^{\pm }=\frac{1}{2}({\widehat{\sigma }}^x\pm \mathrm{i}{\widehat{\sigma }}^y)$, such that the (unnormalized) amplitudes of the ground state in the ${\widehat{S}}^z$ basis follow as

$${\mathsf{\Psi }}_{\mathrm{AKLT}}\left(\mathit{S}\right)=\mathrm{tr}\left({\mathbf{A}}^{S_1}{\mathbf{A}}^{S_2}\cdots {\mathbf{A}}^{S_N}\right).$$

(26)

Since the AKLT point of ${\widehat{H}}_{\mathrm{AKLT}}$ lies in the gapped Haldane phase, $\left|{\mathsf{\Psi }}_{\mathrm{AKLT}}\right.〉$ has finite-ranged magnetic correlations,

$$O_{\ell }^{zz}=\left.〈{\psi }_{\mathrm{AKLT}}\right|{\widehat{S}}_0^z{\widehat{S}}_{\ell }^z\left|{\psi }_{\mathrm{AKLT}}\right.〉\sim e^{-\ell /\xi },\mathrm{with}\xi =\frac{1}{\ln \left(3\right)},$$

(27)

yet it also has an unbroken spin rotation symmetry which is a hallmark of a symmetry protected topological order. Specifically, the string-order parameter

$$O^{\mathrm{string}}=\underset{N\to \infty }{lim}\underset{\ell \to \infty }{lim}\left.〈{\psi }_{\mathrm{AKLT}}\right|{\widehat{S}}_0^z\prod _{j=1}^{\ell -1}{e}^{i\pi {\widehat{S}}_j^z}{\widehat{S}}_{\ell }^z\left|{\psi }_{\mathrm{AKLT}}\right.〉\sim -\frac{4}{9},$$

(28)

reveals the presence of infinite-ranged anti-ferromagnetic correlations. This is evident from the structure of the MPS amplitudes. Any matrix product ${\mathbf{A}}^{\pm }{\mathbf{A}}^0\cdots {\mathbf{A}}^0{\mathbf{A}}^{\pm }=0$, so any configuration containing a ferromagnetic segment, like “+ 0 0 0 +”, with any number of 0’s is not allowed. In contrast, allowed configurations contain only anti-ferromagnetic segments, such as “– 0 + 0 0 0 – + 0”, arising from sequences, like ${\mathbf{A}}^{\pm }{\mathbf{A}}^0\cdots {\mathbf{A}}^0{\mathbf{A}}^{\mp }$, where every ± is partnered with a ∓ separated by with an arbitrary string of 0’s.

Despite its simple MPS representation it is surprisingly non-trivial to capture the non-commutative matrix products making up the AKLT amplitudes with an NQS. Direct conversion to an NQS from the MPS representation gives two reasons why it must contain long-ranged hidden units. First, it has been shown [38] that any short-ranged translationally invariant NQS cast into a MPS form by mapping hidden units into virtual bonds has $\mathbf{A}$ matrices that are at most rank-1. Since the matrix ${\mathbf{A}}^0$ in the AKLT state is rank-2, it fails this condition. Second, if we divide the chain into a sequence of three contiguous parts $a,b,c$, once we make b larger than the longest range of any hidden unit, so no hidden unit connects to visible units in both a and c, then the NQS amplitudes factorize as

$${\psi }_{\mathrm{sr}-\mathrm{NQS}}({\mathit{S}}_a,{\mathit{S}}_b,{\mathit{S}}_c)={\psi }_{ab}({\mathit{S}}_a,{\mathit{S}}_b){\psi }_{bc}({\mathit{S}}_b,{\mathit{S}}_c),$$

(29)

implying that ${\mathit{S}}_a$ and ${\mathit{S}}_c$ are uncorrelated [34]. The AKLT amplitudes ${\mathsf{\Psi }}_{\mathrm{AKLT}}\left(\mathit{S}\right)$ do not satisfy this property since region b can be any length of 0’s, and there will always be non-zero amplitudes “+ 0 0 0 ...0 –” and “–  0 0 0 ...0 +” encoding string order correlations that are not factorizable.

It has been previously found that the AKLT state in the $S^z$ basis requires an NQS with $M\sim O\left(N^2\right)$ long-ranged hidden units [36], and this was borne out in numerical calculations for small systems. In Appendix C, we explicitly construct an NQS for the $S^z$ basis AKLT amplitudes using $M=2N^2+N\lfloor \frac{1}{2}(N-1)\rfloor +1$ hidden units, many of which are extensive over the system. The $O\left(N^2\right)$ scaling can be readily understood as a consequence of having hidden units that each eliminate disallowed configurations, such as “± 0 0 0 ±”, and impose the sign for allowed configurations, such as “± 0 0 0 ∓”, for all N separations and N translations over the system. This is rather less efficient than the compact spin-$\frac{1}{2}$ NQS found for Jastrow, graph and stabilizer states in ref. [42]. The AKLT state can be expressed with $O\left(N\right)$ hidden units but at the expense of needing a 2-layer DBM network [38] that cannot in general be exactly sampled, complicating its use numerically. However, as we saw for the AFH model numerical results, the hidden unit complexity is basis dependent. Surprisingly, we will show next that an exact $M\sim O\left(N\right)$ spin-1 NQS representation of the AKLT state is obtained in the $xyz$ basis.

5.1. Exact Spin-1 NQS for AKLT State in the $xyz$ Basis

The AKLT state provides an instructive example of how a single spin basis change can significantly alter the amplitude structure. Transforming the MPS representation into the $xyz$ basis yields matrices

$$\begin{array}{ccc}\hfill {\mathbf{B}}^x& =& \frac{1}{\sqrt{2}}({\mathbf{A}}^{+}-{\mathbf{A}}^{-1})=\frac{1}{2}{\widehat{\sigma }}^x,{\mathbf{B}}^y=\frac{-\mathrm{i}}{\sqrt{2}}({\mathbf{A}}^{+1}+{\mathbf{A}}^{-1})=\frac{1}{2}{\widehat{\sigma }}^y,{\mathbf{B}}^z=-{\mathbf{A}}^0=\frac{1}{2}{\widehat{\sigma }}^z,\hfill \end{array}$$

and, thus, renders the amplitudes into products of Pauli matrices

$${\mathsf{\Psi }}_{\mathrm{AKLT}}\left(\mathbf{\alpha }\right)=\mathrm{tr}\left({\mathbf{B}}^{{\alpha }_1}{\mathbf{B}}^{{\alpha }_2}\cdots {\mathbf{B}}^{{\alpha }_N}\right),$$

where $\mathbf{\alpha }=({\alpha }_1,\cdots ,{\alpha }_N)$ with ${\alpha }_j=\{x,y,z\}$ label the $xyz$ basis. As expected, there is no change in the complexity/internal dimension of the MPS representation.

The structure of the amplitudes ${\mathsf{\Psi }}_{\mathrm{AKLT}}\left(\mathbf{\alpha }\right)$ in the $xyz$ basis is significantly simpler than ${\mathsf{\Psi }}_{\mathrm{AKLT}}\left(\mathit{S}\right)$ in the $S^z$ basis. Amplitudes are now evaluated by tracking the anticommutations of Pauli matrices required to make the matrices of each the type form a contiguous sequence, e.g., $xx\cdots xyy\cdots yzz\cdots z$, and then reducing the product repeatedly via ${\left({\widehat{\sigma }}^{\alpha }\right)}^2=𝟙$. The resulting matrix trace is non-zero only when the overall product is $\propto 𝟙$, and so all non-zero amplitudes have an equal magnitude. Depending on whether N is even or odd, this condition requires that there is either an even or odd number of $x,y$ and z’s in any configuration string, respectively. Using this, we arrive at the following result:

Theorem 1

(AKLT state $xyz$ NQS). The AKLT state in the $xyz$ spin basis has an exact spin-1 NQS representation requiring $M=2N$ hidden units.

Proof. 

We establish this result using a direct and intuitive construction for ${\mathsf{\Psi }}_{\mathrm{AKLT}}\left(\mathbf{\alpha }\right)$ in which hidden units are devised to implement the nodal structure and sign structure of this state. The rules governing the amplitudes are as follows:

• To implement the parity constraint on the number of $x,y$ and z’s in any configuration string $\mathbf{\alpha }$, we introduce the following $2\times 3$ coupling matrices:

  $${\mathbf{C}}_{xy}=\left[\begin{array}{ccc}\hfill 1& \hfill 1& 1\\ \hfill -1& \hfill -1& 1\end{array}\right],{\mathbf{C}}_{yz}=\left[\begin{array}{ccc}1& \hfill 1& \hfill 1\\ 1& \hfill -1& \hfill -1\end{array}\right].$$
  By defining two hidden units from these matrices as ${{\rm Y}}_{xy}=\{{\mathbf{C}}_{xy},{\mathbf{C}}_{xy},\cdots ,{\mathbf{C}}_{xy}\}$ and ${{\rm Y}}_{yz}=\{{\mathbf{C}}_{yz},{\mathbf{C}}_{yz},\cdots ,{\mathbf{C}}_{yz}\}$, we arrive at the product filter

  $${{\rm Y}}_{xy}\left(\mathbf{\alpha }\right){{\rm Y}}_{yz}\left(\mathbf{\alpha }\right)=\left(1+{(-1)}^{\#x{}^{\prime }\mathrm{s}+\#y{}^{\prime }\mathrm{s}}\right)\left(1+{(-1)}^{\#y{}^{\prime }\mathrm{s}+\#z{}^{\prime }\mathrm{s}}\right),$$

  in which the hidden units cancel out any strings $\mathbf{\alpha }$ that have odd numbers of both x’s and y’s, and y’s and z’s, respectively. Together, these hidden units completely establish for any N the nodal structure of the AKLT state amplitudes in this basis.
• To reproduce the sign structure arising from anticommuting Pauli matrices into a contiguous sequence, we require two types of hidden units. The first type of hidden unit uses a conditional coupling matrix for the local state $\left|x\right.〉$

  $${\mathbf{C}}_x=\left[\begin{array}{ccc}0& \hfill 1& 1\\ 1& \hfill 0& 0\end{array}\right],$$

  along with ${\mathbf{C}}_{yz}$ to define a hidden unit of the form

  $${{\rm Y}}_{x\left(k\right)}=\{{\mathbf{C}}_{yz},\cdots ,{\mathbf{C}}_{yz},\stackrel{k}{\overbrace{{\mathbf{C}}_x}},\mathbf{I},\cdots ,\mathbf{I}\},$$

  where ${\mathbf{C}}_x$ appears in the kth position in the sequence. The action of ${{\rm Y}}_{x\left(k\right)}$ is to induce on a configuration a factor ${(-1)}^{\#y{}^{\prime }\mathrm{s}+\#z{}^{\prime }\mathrm{s}}$ between site $k-1$ and the left boundary, conditional on site k being in state $\left|x\right.〉$. This is the sign that would occur if a ${\sigma }^x$ matrix was anticommuted to this boundary through the corresponding product of Pauli matrices. The second type of hidden unit uses two further coupling matrices

  $${\mathbf{C}}_y=\left[\begin{array}{ccc}1& \hfill 1& 1\\ 1& \hfill -1& 1\end{array}\right],{\mathbf{C}}_z=\left[\begin{array}{ccc}1& 1& 0\\ 0& 0& 1\end{array}\right],$$

  defining a hidden unit of the form

  $${{\rm Y}}_{z\left(k\right)}=\{\mathbf{I},\cdots ,\mathbf{I},\stackrel{k}{\overbrace{{\mathbf{C}}_z}},{\mathbf{C}}_y,\cdots ,{\mathbf{C}}_y\},$$

  where ${\mathbf{C}}_z$ appears in the kth position in the sequence. The action of ${{\rm Y}}_{z\left(k\right)}$ is to induce on a configuration a factor ${(-1)}^{\#y{}^{\prime }\mathrm{s}}$ between site $k+1$ and the right boundary, conditional on site k being in state $\left|z\right.〉$. This is the sign that would occur if a ${\widehat{\sigma }}^z$ matrix was anticommuted to the right boundary through the corresponding product of Pauli matrices, assuming that any ${\widehat{\sigma }}^x$’s have already been anticommuted to the left boundary. To capture all locations k for both types, thus, requires $2(N-1)$ hidden units which entirely establish the sign structure of the AKLT state amplitudes in this basis.

This gives a total of $M=2N$ hidden units. □

The resulting amplitude-wise product decomposition of the AKLT into hidden unit correlators is depicted Figure 7 for $N=6$.

5.2. Analytic Example—AKLT Unary Stabilizer State

An explicit NQS construction for the AKLT state in the $xyz$ basis has been given before in Ref. [40] using unary encoded cells of $\left\{\mathrm{a},\mathrm{b},\mathrm{c}\right\}$ spin-$\frac{1}{2}$’s. Their construction involves initializing the b spin-$\frac{1}{2}$’s in state $\left|+\right.〉$ while entangling the a and c spin-$\frac{1}{2}$’s between adjacent unary cells in the state $\left|\uparrow \right.〉\left|\uparrow \right.〉+\left|\downarrow \right.〉\left|\downarrow \right.〉$. As a tensor network, this is represented as

(30)

where each box is a unary cell, and H denotes the Hadamard matrix. Each unary cell then has the following unitary applied to it

(31)

where S is the phase gate [40] (we have switched the controlled-NOT gates in the circuit given in Figure 5b of ref. [40] into controlled-Z gates here to expose the graph state equivalence). Putting these pieces together, the unary encoded spin-$\frac{1}{2}$ state $\left|{\mathsf{\Psi }}_{\mathrm{unary}}\right.〉$ is a stabilizer state constructed by the circuit

(32)

where the top part (above the dashed line) generates a graph state, and the bottom part applies local Clifford gates. Once unary projection to the spin-1 is applied, $\left|{\mathsf{\Psi }}_{\mathrm{unary}}\right.〉$ generates the AKLT state [40].

Previously, in ref. [42], we showed how any stabilizer state for N spin-$\frac{1}{2}$’s can be readily converted into an NQS with $M\le N-1$. Here, we just summarize the basic process. The first step in this conversion is to use the local Clifford equivalence of stabilizer states to graph states to relocate all the non-diagonal Clifford gates to independent vertices of the graph. This conversion takes the simple chain-like graph state and pattern of Clifford gates from Equation (32) and gives the following for $N=12$ spins:

As required, the resulting graph has diagonal Clifford gates on all vertices, except for a small independent set $\mathcal{I}=\{10,12\}$ highlighted. Notice that the three-site translational invariance of $\left|{\mathsf{\Psi }}_{\mathrm{unary}}\right.〉$, mirrored by the initial chain graph, is still formally present in the transformed graph but is now obscured by its highly connected topology. By forming a vertex cover $\mathcal{C}=\{\mathcal{I},1,2,4,5,7,8\}$, we obtain a NQS [42] with $M=8$ hidden units (despite the more complex graph, this is the same number of hidden units required to describe the initial chain graph state as an NQS):

(33)

More generally, for N spin-$\frac{1}{2}$’s, this procedure generates an NQS with $M=2N/3$ hidden units.

After applying the unary projection and contraction process directly to this spin-$\frac{1}{2}$ NQS, we obtain spin-1 NQS tensor network, whose schematic structure is shown here for $N=4$:

(34)

It is evident that pairs of hidden units possess coordination $4,3$ and 2, while one pair gets projected down to a spin-1 visible bias, giving $M=6$ overall. The same structure applies for general N with pairs of equal coordination spanning $N,N-1,\cdots ,2$, giving a spin-1 NQS for the AKLT state in the $xyz$ basis requiring $M=2N-2$ hidden units. This representation is essentially identical to the one presented in Section 5.1, except that the hidden units implementing the nodal structure (the fully connected pair) also contribute to the sign structure, reducing the total hidden unit count by one pair. This raises an interesting question of whether an even more compressed spin-1 NQS representation of the AKLT state is possible. We finish by examining this using direct numerical optimizations.

5.3. Numerical Example—AKLT in $xyz$ and $S^z$ Bases

Although the AKLT state is translationally invariant, the hidden units encoding the sign structure of this solution are neither individually translationally invariant nor do their translates appear. Consequently, we performed VMC calculations with increasing M using both the spin-1 NQS and unary NQS for $N=\{6,8,10,12\}$ and compared against exact diagonalization. As with the AFH model (which is ${\widehat{H}}_{\mathrm{AKLT}}$ with $\beta =0$), earlier in Section 4.5, we considered both the $S^z$ and $xyz$ basis with their corresponding nodal structure enforced by sampling. We also utilize the same sequential growth scheme as we used in the Heisenberg calculations, again confirming the robustness of our qualitative conclusions by performing reruns of the optimization sequence and presenting the best results here.

As we are performing stochastic optimization, intrinsic sampling noise will limit the accuracy to which any formally exact solution can be found. It is, therefore, crucial to quantitatively characterize when exactness may have been reached numerically. For a gapped Hamiltonian, the average energy of an approximate state E can be related to its infidelity with the ground state using [30]

$$1-\mathcal{O}\le \frac{E-E_0}{\delta }=\frac{ϵ}{\delta },$$

(35)

where $\delta$ is the energy gap of the Hamiltonian, and $ϵ=E-E_0$ is the energy deviation. As the ground state energy $E_0$ of the AKLT Hamiltonian is zero, the energy deviation is simply the sampled energy of the state ${ϵ}_{\mathrm{samp}}$. As shown in Figure 8a, even if the exact analytic spin-1 NQS solution is used, ${ϵ}_{\mathrm{samp}}$ fluctuates when using a finite number of samples typical of an optimization step. For all optimizations presented in this paper, we used the following hyperparameters: number of samples per optimization step $n_{\mathrm{samp}}=8000$, number of optimization steps $n_{\mathrm{step}}\sim 5000$. Typically, the full wavefunction and its fidelity are calculated and checked every 1000 steps to gauge whether the solution has converged or requires further optimization. To account for fluctuations caused by a finite $n_{\mathrm{samp}}$ employed throughout the stochastic optimization, therefore, we use, in Equation (35), $ϵ=\Delta ϵ$, the standard deviation of the sampled energy. As shown in Figure 8b, $\Delta ϵ$ vanishes as $1/\sqrt{n_{\mathrm{samp}}}$, and we estimate an algorithmic fidelity resolution of $\mathcal{R}=\Delta ϵ/\delta \approx 1.2\times {10}^{-5}$ for $N=12$ sites below, in which it may be hard to discriminate an exact solution from an extremely good approximate one.

In Figure 9, we show the smooth decrease in $1-\mathcal{O}$ against M for the $S^z$ basis. While the $N=6$ curve in Figure 9a drops below $\mathcal{R}$, this is not attained for larger N shown in Figure 9b–d within the M’s tested. This indicates that no “exact” NQS solutions have been located. As with the AFH model, the spin-1 NQS and unary NQS achieve a similar accuracy verses M, although the former utilizes less variational parameters.

The analogous results in Figure 10 for the $xyz$ basis display remarkable features in comparison. For each N, a sharp drop in $1-\mathcal{O}$ by over 4 orders of magnitude is observed at $M=N-2$ that consistently pushes the infidelity below $\mathcal{R}$. After reaching this point, $1-\mathcal{O}$ versus M plateaus, and subsequent hidden units have negligible bearing on the accuracy of the wavefunction due to statistical fluctuations originating from the stochastic optimization. These features are consistently produced by both the spin-1 NQS and unary NQS, aside from the largest system size $N=12$ in Figure 10d, indicating that the increasing number of redundant variational parameters in the unary encoding is complicating the optimization. The overall behavior of $1-\mathcal{O}$ observed in this basis is strong evidence of a numerically exact solution with $M=N$ (once the 2 hidden units implementing the nodal structure are included). This is substantially smaller than the analytic solutions introduced and is very suggestive of there being a compact exact spin-1 NQS representation of the AKLT state in the $xyz$ basis.

6. Conclusions

We have introduced the most natural and direct generalization of RBM from spin-$\frac{1}{2}$ to spin-1. This necessitated including a quadratic visible bias and a quadratic visible-hidden interaction in the RBM energy function to ensure trivial product state representation, labeling freedom and gauge equivalence to the tensor network formulation. We demonstrated its use numerically for the spin-1 AFH model in both the $S^z$ and $xyz$ bases, illustrating how the choice of basis can affect the accuracy and hidden unit complexity of an NQS representation. Using our spin-1 NQS, we then re-examined how to represent the AKLT state exactly. In the $S^z$ basis, it is known to require $M\sim O\left(N^2\right)$ hidden units, yet, by changing to the $xyz$ basis, we construct an NQS with $M\sim O\left(N\right)$ hidden units.

Numerical VMC calculations have indicated that, by capturing the nodal structure, either implicitly within the sampling or explicitly through the inclusion of extra hidden units, the optimization can find even more efficient constructions for the sign structure. The resulting spin-1 NQS for the AKLT state in the $xyz$ basis requires $M=N$ hidden units in total making it compact. This example raises the interesting possibility of improving the efficiency and accuracy of NQS calculations by including single-spin basis transformations to lower the hidden unit complexity.

Several important open questions remain about NQS representations. In particular, it would be instructive to build representations of classes of bosonic states using multinomial RBMs. In this case, a local Fock basis is typically employed; however, our findings suggest that it could be useful to explore a local basis that breaks the particle number symmetry when describing condensates. Moreover, the elevation of visible units from binary to multinomial raises the question of whether also using multinomial hidden units can enhance the expressiveness of NQS. This has been explored in the context of binary visible units in ref. [54] in an analytical context to precisely represent certain two- and three-body interactions. The use of multinomial hidden units for numerical VMC calculations has been largely unexplored and is the subject of forthcoming work [55] for the Bose Hubbard model in two dimensions.