Approximating power of machine-learning ansatz for quantum many-body states

TL;DR

This paper analyzes the effectiveness of restricted Boltzmann machine neural networks as variational ansatzes for quantum many-body states, revealing their connection to perturbation theory and their performance limitations related to entanglement.

Contribution

It uncovers the relation between RBM structure and perturbation series, improves algorithms for better convergence, and compares RBM with other variational wave-functions.

Findings

RBM ansatz achieves high precision in simple models due to perturbation series connection.

Performance of RBM depends on entanglement properties of target states.

Introducing alternative wave-functions yields comparable results to RBM.

Abstract

An artificial neural network (ANN) with the restricted Boltzmann machine (RBM) architecture was recently proposed as a versatile variational quantum many-body wave function. In this work we provide physical insights into the performance of this ansatz. We uncover the connection between the structure of RBM and perturbation series, which explains the excellent precision achieved by RBM ansazt in certain simple models, demonstrated in the literature. Based on this relation, we improve the numerical algorithm to achieve better performance of RBM in cases where local minima complicate the convergence to the global one. We introduce other classes of variational wave-functions, which are also capable of reproducing the perturbative structure, and show that their performance is comparable to that of RBM. Furthermore, we study the performance of a few-layer RBM for approximating ground states of random, translationally-invariant models in 1d, as well as random matrix-product states (MPS). We find that the error in approximating such states exhibits a broad distribution, and is largely determined by the entanglement properties of the targeted state.

Figures

Equations

$${\Psi }_{RBM}(\{{\sigma }_j\})=\sum _{\{s_i=\pm 1\}}{e}^{{\sum }_j{a}_j{\sigma }_j^z+{\sum }_i{b}_i{s}_i+{\sum }_{ij}{W}_{ij}{s}_i{\sigma }_j^z},$$

$${\Psi }_{RBM}(\{{\sigma }_j\})\propto$$

$$\prod _{i=1}^M\cosh \left(b_i+\sum _j{W}_{ij}{\sigma }_j^z\right).$$

$$H_{TFI}=-\sum _i{\sigma }_i^x-J\sum _i{\sigma }_i^z{\sigma }_{i+1}^z,$$

$$|\psi ⟩=\frac{e^{\hat{W}}|{\psi }_0⟩}{⟨{\psi }_0|e^{{\hat{W}}^{†}}{e}^{\hat{W}}|{\psi }_0{⟩}^{\frac{1}{2}}},$$

$$(H{e}^{\hat{W}}-e^{\hat{W}}E)|{\psi }_0⟩=0.$$

$${\Psi }_{RBM}(\{{\sigma }_j\})=e^{\tilde{W}}|{\tilde{\psi }}_0⟩,$$

$$\Psi_{quad}(\{\sigma_{j}\})=\prod^{M}_{i=1}\left(1+\alpha_{i}\Big{[}b_{i}+\sum_{j}W_{ij}\sigma^{z}_{j}\Big{]}^{2}\right),$$

$${\Psi }_{fix}(\{{\sigma }_j\})=\prod _{i=1}^M\sum _{k=0}^3{\left(b_i+\sum _j{W}_{ij}{\sigma }_j^z\right)}^{2k}.$$

$$H_{ni}=-\sum _i{\sigma }_i^x-h\sum _i{\sigma }_i^z-J\sum _i{\sigma }_i^z{\sigma }_{i+1}^z,$$

$$H_{rand}=\sum _i\sum _{\alpha ,\beta =\{0,x,y,z\}}{J}_{\alpha \beta }{\sigma }_i^{\alpha }{\sigma }_{i+1}^{\beta },$$

$${\Psi }_{RBM}(\{{\sigma }_j\})=\prod _i{f}_{{\sigma }_i..{\sigma }_{i+d-1}}^{(i)},$$

$$A_{j_1..j_{d-1}{k}_1..k_{d-1}}^{{\sigma }_i}=f_{{\sigma }_i{k}_1..k_{d-1}}^{(i)}{\delta }_{{\sigma }_i{j}_1}{\delta }_{k_1{j}_2}\cdot ..\cdot {\delta }_{k_{d-2}{j}_{d-1}},$$

$${\Psi }_{RBM}=\sum _{\{{\sigma }_j\}}Tr(A^{{\sigma }_1}..A^{{\sigma }_N})|\{{\sigma }_j\}⟩$$

$$(H{e}^{\hat{W}}-e^{\hat{W}}E)|{\psi }_0⟩=0.$$

$$\left([H_0,\widehat{W_1}]+V\right)|{\psi }_0⟩=0,$$

$$\left([H_0,\widehat{W_2}]+\frac{1}{2}[[H_0,\widehat{W_1}],\widehat{W_1}]\right)|{\psi }_0⟩=E_2|{\psi }_0⟩,$$

$$\left([H_0,\widehat{W_3}]+[[H_0,\widehat{W_1}],\widehat{W_2}]+\frac{1}{6}[[[H_0,\widehat{W_1}],\widehat{W_1}],\widehat{W_1}]\right)|{\psi }_0⟩=0.$$

$$\frac{1}{2}[[H_0,\widehat{W_1}],\widehat{W_1}]=-4{\lambda }^2\sum _i({\sigma }_i^x+{\sigma }_i^z{\sigma }_{i+1}^x{\sigma }_{i+2}^z).$$

$$[[H_0,\widehat{W_1}],\widehat{W_2}]+\frac{1}{6}[[[H_0,\widehat{W_1}],\widehat{W_1}],\widehat{W_1}]=\frac{16{\lambda }^{3}i}{3}\sum _i\left({\sigma }_i^y{\sigma }_{i+1}^z+{\sigma }_i^z{\sigma }_{i+1}^y\right)-$$

$$4{\lambda }^3\sum _i\left({\sigma }_i^x{\sigma }_{i+1}^z{\sigma }_{i+2}^z+{\sigma }_i^z{\sigma }_{i+1}^z{\sigma }_{i+2}^y+{\sigma }_i^z{\sigma }_{i+1}^x{\sigma }_{i+3}^z+{\sigma }_i^z{\sigma }_{i+2}^x{\sigma }_{i+3}^z\right).$$

$$|\psi ⟩=e^{\hat{W}}|{\psi }_0⟩+O({\lambda }^4),$$

$$\hat{W}=\sum _i\left((\lambda -\frac{2{\lambda }^3}{3}){\sigma }_i^z{\sigma }_{i+1}^z+{\lambda }^2{\sigma }_i^z{\sigma }_{i+2}^z+2{\lambda }^3{\sigma }_i^z{\sigma }_{i+3}^z\right).$$

$${\Psi }_{RBM}(\{{\sigma }_j\})={\left(\prod _{i=0}^3\cosh {\tilde{W}}_i\right)}^N\prod _j(1+w_0{w}_1{\sigma }_j^z{\sigma }_{j+1}^z+w_0{w}_2{\sigma }_j^z{\sigma }_{j+2}^z+w_0{w}_3{\sigma }_j^z{\sigma }_{j+3}^z+w_1{w}_2{\sigma }_{j+1}^z{\sigma }_{j+2}^z+$$

$$+w_1{w}_3{\sigma }_{j+1}^z{\sigma }_{j+3}^z+w_2{w}_3{\sigma }_{j+2}^z{\sigma }_{j+3}^z+w_0{w}_1{w}_2{w}_3{\sigma }_j^z{\sigma }_{j+1}^z{\sigma }_{j+2}^z{\sigma }_{j+3}^z),$$

Full text

Introduction.— Variational methods play an invaluable role in quantum many-body physics, because they allow one to represent exponentially many amplitudes of a many-body wave function using a small number of variational parameters. The choice of the variational ansatz is often motivated by the underlying physics of the system of interest, notable examples being product states, BCS wave function (Bardeen et al., 1957), and Laughlin states (Laughlin, 1983). Broad classes of variational wave functions, such as tensor networks (Orús, 2014; Perez-Garcia et al., 2007; Verstraete et al., 2008), which include matrix product states (MPS) (Perez-Garcia et al., 2007), rely on the low amount of quantum entanglement in ground states of physical systems.

Even though tensor network methods proved remarkably successful (White, 1992; Schollwöck, 2011; Vidal, 2003; Verstraete et al., 2008; Vidal, 2008), it is important to investigate other classes of variational wave functions, including the ones that can capture states with higher entanglement. Recent proposals (Carleo and Troyer, 2017; Saito, 2017a; Gao and Duan, 2017; Nomura et al., 2017; Glasser et al., 2018; Cai and Liu, 2018a; Kaubruegger et al., 2018; Liang et al., 2018; Pastori et al., 2018; Saito, 2017b; Cai and Liu, 2018b) for such ansätze considered variational functions inspired by machine learning (ML). More generally, ML is finding an increasing number of diverse applications in physics, including detection of phase transitions (Ohtsuki and Ohtsuki, 2016; van Nieuwenburg et al., 2017; Carrasquilla and Melko, 2017; Broecker et al., 2017; Ch’ng et al., 2017; Ohtsuki and Ohtsuki, 2017; Schindler et al., 2017; Tanaka and Tomiya, 2017; Venderley et al., 2018; Hu et al., 2017), extraction of relevant degrees of freedom (Torlai and Melko, 2016; Koch-Janusz and Ringel, 2018), and improvement of existing techniques (Huang and Wang, 2017; Liu et al., 2017; Torlai et al., 2018). In this paper we consider the variational ansatz inspired by the restricted Boltzmann machine (RBM) (Carleo and Troyer, 2017). The architecture of RBM – a neural network with a wide range of applications outside quantum physics (Hinton and Salakhutdinov, 2006; Salakhutdinov et al., 2007; Larochelle and Bengio, 2008) – is illustrated in Fig. 1.

The RBM ansatz for the system of $N$ physical spins $\{\sigma_{i}^{z}\}=\pm 1$ is constructed by introducing $M$ hidden spins $\{s_{i}\}=\pm 1$, which gives $\alpha=M/N$ hidden spins per each physical spin. Then, the amplitudes of the variational wave function in the $\{\sigma^{z}_{i}\}$ eigenbasis are obtained by summing over the states of hidden degrees of freedom:

[TABLE]

This wave function is optimised over the variational parameters $\{a_{j}\}$, $\{b_{i}\}$, $\{W_{ij}\}$ to yield the lowest-energy wave function for a given Hamiltonian.

The summation over the hidden-spin states in Eq. (1) can be performed explicitly, giving (up to a normalization factor):

[TABLE]

This property facilitates Monte Carlo sampling, since any quantum amplitude can be computed by evaluating the function (2).

The utility of the RBM ansatz (2) is being actively explored (Deng et al., 2017a, b; Chen et al., 2018; Jia et al., 2019; Zhang et al., 2018; Lu et al., 2018). Among its advantages are applicability in any number of spatial dimensions, and the ability to describe states with high degree of entanglement. It was found to approximate the ground states (GS) of certain Hamiltonians (1D Ising, 1D and 2D Heisenberg) (Carleo and Troyer, 2017) with a remarkably high precision, given relatively few variational parameters. Moreover, it can exactly represent certain topological states (1D cluster state, 2D and 3D toric code) (Deng et al., 2017b).

An important question regarding the representation power of tensor-network states by RBM ansatz is intimately related to the entanglement properties of the latter. The entanglement entropy of a general RBM state obeys volume-law scaling (Deng et al., 2017a), but finite-range RBM states (that is, states in which a hidden spin can have non-zero interaction $W_{ij}$ only with $d$ contiguous physical spins) obey area-law in any number of dimensions (Deng et al., 2017a), and can be represented by an MPS with a finite bond dimension $D$ [Chen et al., 2018]. However, some MPS states, AKLT state being an example, cannot be approximated by a finite-range RBM (Chen et al., 2018). The infinite-range RBM can arbitrary well approximate any MPS, but the number of required hidden units is exponential in the bond dimension of MPS (Huang and Moore, 2017), which is impractical for numerical applications. Thus, it is important to understand how well shallow RBM states, with relatively low density of hidden spins $\alpha$, which can be efficiently optimized, can approximate generic MPS states.

Here, we provide several results regarding the performance, representation power, and the uniqueness of the RBM ansatz. First, by connecting RBM to perturbation theory, we explain the surprisingly good performance of RBM with few variational parameters for certain models, reported in Ref. (Carleo and Troyer, 2017). To that end, we show that RBM with $\alpha=1$ already captures several orders of the perturbation series for those models. The connection to the perturbation theory naturally suggests an improvement for the optimisation algorithm, which we explore. Second, this observation points to a whole class of variational wave functions, similar to RBM; we demonstrate that their performance in simple models is comparable to that of RBM. Third, we investigate the performance of RBM for general local Hamiltonians, picked at random, as well as its ability to approximate random MPSs in the practically interesting case of $\alpha=1$. By applying the ansatz to various random realisations of these systems, we show that the performance of RBM is universally determined by the entanglement entropy and range of correlations.

RBM and perturbation theory.— To illustrate the connection of the structure of RBM to perturbation series we first focus on a solvable, 1D transverse-field Ising (TFI) model:

[TABLE]

where $\sigma_{i}^{\beta},\beta=x,y,z$ are the standard Pauli operators.

We consider two points in the phase diagram: $J=0.5$ (paramagnetic phase), and $J=2$ (ferromagnetic phase). The GS wave function at these points can be presented in the following form, using perturbation theory:

[TABLE]

where $|\psi_{0}\rangle$ is the unperturbed wave function: polarized along $x$-axis or $z$-axis for the paramagnetic and ferromagnetic case, respectively. The term $\hat{W}=\sum_{n}\hat{W}_{n}$ has a perturbative structure, with the $n$-th term accounting for the $n$-th order correction in a perturbation series. The terms $\hat{W}_{n}$ become less and less local as $n$ is increased. Its explicit form is obtained by iteratively solving the Schrödinger equation

[TABLE]

The RBM wave function (2) can be rewritten in a form similar to Eq. (4):

[TABLE]

where $\tilde{W}=\sum_{j}a_{j}\sigma^{z}_{j}+\sum^{M}_{i=1}\ln\cosh\left(b_{i}+\sum_{j}W_{ij}\sigma^{z}_{j}\right)$ and $|\tilde{\psi}_{0}\rangle$ stands for the state polarised along the $x$-axis.

We note that for the paramagnetic case $\hat{W}$ can be expressed solely in terms of $\{\sigma^{z}_{i}\}$ operators and the [math]th order wave function $|\psi_{0}\rangle$ coincides with $|\tilde{\psi}_{0}\rangle$ (see Appendix for more details). This allows us to build a map between the perturbation series and the RBM ansatz at the paramagnetic point. By choosing $\alpha=1$, $a_{j}=0$, $b_{i}=0$ and $W_{ij}=\sum_{j^{\prime}=0}^{3}\delta_{i,j+j^{\prime}}\text{arctanh}\ w_{j^{\prime}}$ with $w_{0}=1-J^{4}$, $w_{1}=J-2J^{3}$, $w_{2}=J^{2}+J^{3}$, $w_{3}=2J^{3}$ we are able to reproduce the first three orders of the perturbation series exactly.

Clearly, the RBM ansatz is not capable of capturing all orders of perturbation theory. The $n$-th order term in the perturbative expansion, $\hat{W}_{n}$, generally contains operators with the support that scales as $O(n)$. Thus, the number of all possible terms at $n$th order scales exponentially with $n$, and exponentially many variational parameters would be needed to exactly reproduce that. However, even a finite number of terms is sufficient to approximate GS energy with a high precision. This is illustrated in Fig. 2. At low orders the curve describing RBM precision of the paramagnetic GS energy approximation (blue) follows closely the expansion in powers of $J$ of the exact energy (purple), computed using Jordan-Wigner transformation. At high orders ($\geq 8$), the RBM performance becomes worse, which suggests that higher order corrections are captured only approximately.

In Fig. 2 we relate a given RBM ansatz to a certain order of perturbation theory. We achieve this by taking the range $d$ of interactions $W_{ij}$ to be equal to the largest support among the irreducible terms that contribute to the corresponding order of the perturbative expansion of the energy. For example, the second order correction to the energy is determined by the first order correction to the GS wave function in the case of TFI model, which gives $d=2$ ($d=1$) for the paramagnetic (ferromagnetic) point.

The explicit construction used above does not directly carry over to the ferromagnetic case, since in the unperturbed GS $|\psi_{0}\rangle$ spins are aligned along $z$-axis rather than $x$-axis and consequently $\hat{W}$ contains the spin-flip terms ($\sigma^{x}_{i}$ or $\sigma^{y}_{i}$), absent in the RBM wave function (2). A straightforward way to circumvent this is to rotate the basis to make the spins in the unperturbed GS polarized along the $x$-axis. In particular, this allows us to get a relative energy error $\sim\!\!10^{-6}$ for the ferromagnetic point of 1D TFI, which is one order of magnitude smaller than the result obtained without rotation of the basis Carleo and Troyer (2017) (yellow dashed line in Fig. 2).

Practically, the parametric argument suggests that in RBM states with relatively low $\alpha\gtrsim 1$, there should be enough variational parameters to capture several lowest orders of the perturbative expansion. There is however a tradeoff: RBM states with larger values of $\alpha$, while forming a larger variational manifold, are harder to optimize efficiently, as one can get stuck in local minima. To enforce perturbative structure we slightly modify the original optimisation algorithm (Carleo and Troyer, 2017) in the following way: we gradually increase the range of variational couplings $W_{ij}$, at each step starting with the optimal parameters obtained in the previous step. We found that this modified procedure has a faster convergence, yielding lower error that the original procedure already at range $d=3$.

A family of ansätze.— Another interesting question concerns the uniqueness of the RBM ansatz. The above discussion suggests that it belongs to a whole class of variational wave functions which can reproduce perturbative expansions. To demonstrate this, we investigate two examples below. In the first example, each hidden spin generates a quadratic polynomial:

[TABLE]

and coefficients $\alpha_{i}$ are variational parameters. The second trial ansatz is sixth-order polynomial with only even powers and fixed coefficients

[TABLE]

We restrict ourselves to even powers only, since we found such ansätze to yield a better performance.

The approximative power of these ansätze is summarized in the Fig. 3, where we plot the accuracy of the GS energy approximation as a function of the hidden-unit density, either using the original algorithm or the modified one. The numerical results are presented for the system of $N=18$ spins and four Hamiltonians, namely, paramagnetic ($J=0.5$), ferromagnetic ($J=2$), critical ($J=1$) and a non-integrable one, defined as

[TABLE]

where we choose $J=1$ and $h=0.5$. All three ansätze exhibit similar numerical performance, which indicates that the RBM ansatz belongs to a wider class of functions that are capable of capturing local correlations.

General Hamiltonians.— While TFI model provides a good test for the variational wave functions, it is important to test RBM for more general Hamiltonians. As an example, we consider a family of translationally invariant Hamiltonians with random nearest neighbour spin-spin interactions:

[TABLE]

with $J_{\alpha\beta}$ being independent random coefficients uniformly distributed in the range $[-0.5,0.5]$. We post-select Hamiltonians with translationally invariant GS wave functions and use the RBM ansatz with $\alpha=1$ and translationally invariant variational parameters to approximate the GS energy of these Hamiltonians.

By testing the RBM performance for different realizations of the Hamiltonian (10), we discovered that the error in approximating the GS energy has a very broad distribution. To analyze the nature of this spread, we also studied the properties of the exact GS of the corresponding Hamiltonians (computed using exact diagonalization of system with $N=18$). Namely, we study how the RBM performance depends on the entanglement entropy and the correlation length 111If $K_{i}(n)=\langle\sigma_{1}^{i}\sigma_{1+n}^{i}\rangle_{c}$, then the correlation length $l$ is defined by $l=\textrm{max}_{i=\{x,y,z\}}\sum_{n}\textrm{dist}(1,1+n)|K_{i}(n)|/\sum_{n}|K_{i}(n)|$, where $\textrm{dist}(i,j)$ is the distance between to spins at sites $i,j$ with the periodically boundary conditions being taken into account (Fig. 4). The positive correlation of the error with both quantities is evident, indicating that RBM performs best for states with short-ranged correlations and low entanglement, as expected.

RBM and MPS.— We can further understand the interplay of RBM and entanglement by studying the relation of the ansatz to MPSs. Let us first provide an exact and transparent mapping from finite-range RBM to MPS in 1D. The RBM ansatz with coupling range $d$ can be written as

[TABLE]

where $f^{(i)}_{\sigma_{i}..\sigma_{i+d-1}}$ is a function acting on $d$ neighbouring spins. For example, to write the RBM wave function with $\alpha=1$ in this form, we can choose $f^{(i)}_{\sigma_{i}..\sigma_{i+d-1}}=e^{a_{i}\sigma_{i}^{z}}\cosh(b_{i}+\sum^{d-1}_{j=0}W_{i,i+j}\sigma^{z}_{i+j})$. Next, we construct a mapping from this form to an MPS form. We can associate the tensor $A^{\sigma_{i}}_{i-d,..,i+d-1}$ with each function $f^{(i)}_{\sigma_{i}..\sigma_{i+d-1}}$ using the following rule

[TABLE]

where the first $d-1$ indices ($j_{1}..j_{d-1}$) of the tensor should be considered as one index taking $2^{d-1}$ values, which should be contracted with the $(i-1)$th tensor and the last $d-1$ indices ($k_{1}..k_{d-1}$) are contracted with the $(i+1)$th tensor. Then the RBM-wave function can be written in the form

[TABLE]

Our argument makes use of the fact that the spin system either has periodic boundary conditions or is infinite, however, this argument can be effortlessly generalized to finite size systems with boundaries. Also the slightly modified argument relates RBM to MPS in higher dimensions.

In contrast to Ref. [Chen et al., 2018], where the mapping was made to an MPS with a bond dimension $D$ that is exponential in the number of pairs of the connected hidden unit and the physical spin that have a given spin in between ($\sim\leavevmode\nobreak\ \!\!e^{d^{2}}$), we get an exponential ($D\sim e^{d}$) scaling with the range of RBM.

To study the inverse mapping, we use RBM with $\alpha=1$ to approximate MPSs with a given bond dimension. We generate random translationally invariant MPSs $|\Psi_{MPS}\rangle$ with $D=2,4,8$ for a system of $N=18$ spins 222We construct MPSs as products of identical random matrices that have independent random entries with both real and imaginary parts being uniformly distributed in the interval $[-0.5,0.5]$.. Then we minimize the quantity $1-|\langle\Psi_{MPS}|\Psi_{RBM}\rangle|^{2}$, which is plotted against the entanglement entropy and the correlation length in Fig. 4. One can see a clear dependence of the error on the entanglement entropy, which is in agreement with the results obtained for random Hamiltonians above.

Conclusions.— In this paper we considered the RBM “black box” from the physical prospective. We pointed out an explicit connection between the RBM ansatz and perturbative description of ground states of gapped models, which explains the remarkable accuracy of few-layer RBM for simple models such as the TFI model. [Carleo and Troyer, 2017]. In some cases, several orders of perturbation series can be exactly captured by the RBM wave function. Even when this is not the case, the mere existence of the perturbative series allowed us to introduce a simple modification of the optimization algorithm of RBM with an improved performance. Furthermore, the intuitive connection with the perturbation theory helped us to introduce a whole family of ansätze that have a performance similar to that of RBM. This suggests that RBM is not a unique ansatz, but rather belongs to a broad class of similarly powerful variational wave functions.

Finally, we have investigated the performance of practical, few-layer RBM states, in approximating ground states of random (but translationally invariant) local spin Hamiltonians. We found that the approximation error has a broad distribution, and the success of few-layer RBM is mostly determined by the degree of entanglement in the ground state. It is natural to expect that states with low entanglement can be well approximated using just few orders of perturbation theory, which RBM can capture.

We leave for the future work the detailed investigation of how well infinite-range RBM can reproduce critical many-body states. Another open question is whether the volume-law entanglement that RBM states generally have can be useful for approximating, e.g., excited states and entanglement spreading in many-body systems.

Acknowledgement.— This work was supported by the Swiss National Science Foundation. We thank G. Carleo for helpful correspondence. The numerical simulations were performed on the HPC cluster Baobab.

Appendix A Appendix: An exact mapping of the third-order perturbative wave function to the RBM ansatz

As it is stated in the main text (Eq. (4)) $\hat{W}$ is determined by

[TABLE]

Let us write $\hat{W}$ and the energy $E$ in the form $\hat{W}=\hat{W}_{1}+\hat{W}_{2}+\hat{W}_{3}+O(J^{4})$ and $E=E_{0}+E_{1}+E_{2}+E_{3}+O(J^{4})$ to outline their perturbative structure.

To simplify analysis we want to fix the freedom in choice of $\hat{W}$. It could be chosen to depend solely on $\{\sigma^{z}_{i}\}$. This is the case, since an action of any operator on $|\psi_{0}\rangle$ that is the product of $\{\sigma^{\alpha}_{i}\}$ results only in spin flips at certain positions. Thus, to cancel such terms we need to put $\sigma^{z}_{i}$ in the corresponding positions of $\hat{W}$ . When commuted with $H_{0}$, some $\sigma^{z}_{i}$ would be changed to $\sigma^{y}_{i}$, but the resulting action on $|\psi_{0}\rangle$ stays the same up to a phase ($\sigma^{z}_{i}|\psi_{0}\rangle=i\sigma^{y}_{i}|\psi_{0}\rangle$).

Two important simplifications follow from this observation . Namely that $[V,\hat{W}]=0$ at any order and $e^{\hat{W}_{n}+\hat{W}_{m}}=e^{\hat{W}_{n}}e^{\hat{W}_{m}}$. Also, a simple property of the TFI Hamiltonian is that all odd energy corrections are zero. Having said these the equations that determine the terms of interest in $\hat{W}$ can be written as

[TABLE]

To have more compact notation we denote $\lambda\equiv\frac{J}{4}$. If we take $\hat{W}_{1}=\lambda\sum_{i}\sigma^{z}_{i}\sigma^{z}_{i+1}$, then $[H_{0},\hat{W_{1}}]=2i\lambda(\sum_{i}\sigma^{y}_{i}\sigma^{z}_{i+1}+\sum_{i}\sigma^{z}_{i}\sigma^{y}_{i+1})$. Therefore,the term that should be cancelled in the second order is

[TABLE]

We choose $\hat{W}_{2}=\lambda^{2}\sum_{i}\sigma^{z}_{i}\sigma^{z}_{i+2}$ and get an extensive energy correction $E_{2}=-4\lambda^{2}\sum_{i}1$, then $[H_{0},\hat{W_{2}}]=2i\lambda^{2}(\sum_{i}\sigma^{y}_{i}\sigma^{z}_{i+2}+\sum_{i}\sigma^{z}_{i}\sigma^{y}_{i+2})$ . The third order interaction is then given by

[TABLE]

This gives $\hat{W}_{3}=2\lambda^{3}\sum_{i}(\sigma^{z}_{i}\sigma^{z}_{i+3}+\sigma^{z}_{i}\sigma^{z}_{i+1}/3)$ and the total wave function up to normalization could be written as

[TABLE]

with

[TABLE]

At this point is is clear how RBM ansatz could be constructed in order to recover perturbative expansion, i.e. Taylor expansion of each cosine should reproduce local terms in the exponent of 21. Since we need terms only even in operators $\sigma^{z}$, it is natural to take parameters $a_{i}$, $b_{j}$ of RBM ansatz to zero. The parameters $W_{ij}$ are chosen to be translationally invariant, i.e. dependent only on difference between the position of a hidden-unit and the physical spin $W_{ij}=\tilde{W}_{i-j}$. Since the largest support of the terms in $\hat{W}$ is four, we restrict ourselves to four non-zero parameters $\tilde{W}_{i}$, $\{i=0,..,3\}$. Applying all the assumptions to the RBM wave function (Eq.(2)), we may write ansatz in the form

[TABLE]

where $w_{i}=\tanh\tilde{W_{i}}$. Since we assume this expression to have the perturbative structure , the terms containing Pauli matrices should be small. Thus, we can take approximate logarithm of this expression to get the exponential form as of Eq. 21 up to normalization. Then we could easily identify the value of parameters $w_{i}$ that reproduce the perturbative expansion to the third order: $w_{0}=1$ (it is possible with an exponential precision for finite $\tilde{W}_{0}$ ), $w_{1}=\lambda-2\lambda^{3}$, $w_{2}=\lambda^{2}+\lambda^{3}$ and $w_{3}=2\lambda^{3}$.