Backflow Transformations via Neural Networks for Quantum Many-Body Wave-Functions
Di Luo Bryan K. Clark




















 






 

 









    *#1*/Institute for Condensed Matter Theory and Department of Physics, University of Illinois at Urbana-Champaign, IL 61801, USAInstitute for Condensed Matter Theory and Department of Physics, University of Illinois at Urbana-Champaign, IL 61801, USA
    
    Obtaining an accurate ground state wave function is one of the great challenges in the quantum many-body problem.  In this paper, we propose a new class of wave functions, neural network backflow (NNB).  The backflow approach, pioneered originally by Feynman <cit.>, adds correlation to a mean-field ground state by transforming the single-particle orbitals in a configuration-dependent way.  NNB uses a feed-forward neural network to learn the optimal transformation via variational Monte Carlo.  NNB directly dresses a mean-field state, can be systematically improved and directly alters the sign structure of the wave-function.  It generalizes the standard backflow<cit.> which we show how to explicitly represent as a NNB.  We benchmark the NNB on Hubbard models at intermediate doping finding that it significantly decreases the relative error, restores the symmetry of both observables and single-particle orbitals, and decreases the double-occupancy density. Finally, we illustrate interesting patterns in the weights and bias of the optimized neural network.
    
    Backflow Transformations via Neural Networks 
     for Quantum Many-Body Wave-Functions
        Bryan K.  Clark
        April 14, 2020
    ===================================================================================
    
    Introduction.—A key question in strongly correlated quantum systems is to obtain an approximation for the ground state wave function. This is especially important for Fermion systems in two or more dimensions where only approximate or exponentially costly methods for evaluating observables of quantum systems exist.  Early attempts for writing down variational Fermion wave-functions, such as Slater determinants <cit.> and BCS wave-functions<cit.>, focused on finding the ground state of a mean field Hamiltonian which best matched the interacting ground state.  Since these early attempts more sophisticated wave-functions have been developed which dress these mean-field starting points including Slater-Jastrow <cit.>, Slater-Jastrow-Backflow <cit.> and iterative backflow <cit.> which has recently been described as a non-linear network <cit.>.  These wave-functions have the advantage that the mean-field starting point can directly incorporate the basic physics of the problem.  
    
    Instead of starting from a dressed mean-field, many other classes of wave-functions are parameterized by a tuning parameter D which interpolates from a trivial state at small D to a universal wave-function spanning the entire Hilbert space at exponential D. Examples of such wave-functions include matrix-product states<cit.>,  other forms of tensor networks<cit.>, Huse-Elser states<cit.>, and string-bond states <cit.>.
    Recently, wave-functions based on neural network primitives, such as restricted Boltzmann machines (RBM) and feed forward neural network(FNN), have been introduced with similar universal properties<cit.>:  as the number of hidden neurons increases, the neural network state can represent all probability distribution although may require complex weights to represent the sign structure of Fermion wave-functions.  A recent attempt to incorporate RBM into Fermion states by using the RBM as a more general Jastrow <cit.> shows promise but was still restricted to the sign-structure of the underlying mean-field ansatz. Even though general neural networks could alter the sign-structure, it may struggle with capturing the underlying mean-field physics both in terms of the number of neurons required as well as optimization. 
    
    
    
    In this work, we propose a new class of wave-functions, the Neural Network Backflow (NNB), which dresses a  mean-field wave-function, can make changes to the sign structure directly, and can be systematically improved by increasing the number of hidden neurons and is able to be made theoretically exact in the limit of enough neurons.  To accomplish this we use a feed-forward neural network (FNN), not in the standard approach of returning a wave-function amplitude, but instead to transform the single particle orbitals in a configuration dependent way; these orbitals are then used in the mean-field wave-function.  Wave-functions with configuration-dependent orbitals are known as a backflow wave-function<cit.>. 
    
    
    Background.—Mean field theory - approximating the ground state of a quartic Hamiltonian by the ground state, ψ_MF, of a quadratic Hamiltonian - is a powerful first step to understanding correlated quantum systems.    Various mean fields lead to  different types of ground states including Slater determinants, 
    
        ψ_SD(r)  =     [M^SD,↑] [M^SD,↓]; 
        
        M_ik^SD,σ =    ϕ_kσ(r_iσ)
    
    and Bogoliubov de Gennes wave-functions, 
    
        ψ_BDG(r)=    [Φ] 
        Φ_ij =    ∑^N_k,l=1ϕ_k↑(r_i,↑)S_klϕ_l↓(r_j,↓)
    
    where ϕ_kσ is the k'th single particle orbital and r_i,σ is the position of the i'th particle of spin σ. Eq. (<ref>) only takes the occupied orbitals while Eq. (<ref>) is summed over both occupied and unoccupied orbitals.
    
    Mean field states are uncorrelated by construction.  The simplest way to capture correlation physics is through the introduction of a Jastrow giving ψ_Jastrow(r) =  exp[-U(r)]ψ_MF(r), where U(r) is an arbitrary function. In this work, we always use a charge Jastrow U(r)=1/2∑_i,j v_ij n_i n_j, where n_i is the charge density, v_ij is the variational parameters. While Jastrow factors can introduce many-body correlations, they can't modify the mean-field's sign structure. One approach to add additional sign-structure modifying correlation  is through a backflow correction <cit.>  which introduces correlations by having the single-particle orbitals act on a configuration-dependent quasi-particle position.  
    On the lattice, the backflow approach instead uses a configuration-dependent mean-field <cit.> - i.e. the quadratic Hamiltonian or single-particle orbitals  ϕ_k,σ^b(r_i;r) depend not only on the position r_i but on all other electron positions r; our NNB builds on top of this formulation of backflow.
    
    Neural Network Backflow.—The NNB uses a FNN to modify the single particle orbitals for a spin σ,
    
        ϕ_kσ^b(r_i,σ;r)= ϕ_kσ(r_i,σ) +  a^NN_ki,σ(r)
    
    where each value of  a^NN_ij,σ is represented by an output neuron of the FNN.  We use one neural net for each of σ∈{↑, ↓}. This is to be contrasted with the standard backflow<cit.> parameterization, 
    
        ϕ_kσ^b(r_i,σ;r)     =    ϕ_kσ + ∑_jη_ij,σϕ_kσ(r_j,σ) 
        η_ij,σ     =     t D_i H_j θ_|i-j|,σ
    
    with D_i = n_i,↑n_i,↓, H_i=(1-n_i,↑)(1-n_i,↓). θ_1,σ and θ_2,σ are the only non-zero variational parameters. 
    
    
    Interestingly, the backflow transformation of Eq. (<ref>), can be represented as a neural network for a^NN_ij,σ(r) with three hidden layers and a linear number of neurons; an explicit construction will be given in the next section. This ensures that there exists a three layer neural network which is at least as good as the standard backflow transformation.
    
    
    We consider two NNB wave-functions, Ψ_SN and Ψ_PN, implemented on top of a Slater Determinant and BCS pairing wave-functions respectively.  The neural nets used in these wave-functions are similar although Ψ_SN has outputs which only correspond to the occupied orbitals, while the outputs of Ψ_PN correspond to all the orbitals. In addition, for Ψ_SN there are only two neural nets (one for each of the spin-up and spin-down orbitals) while for Ψ_PN there is an additional neural net used to generate a system dependent S_kl(r). This is implemented by letting S_kl(r)=S_kl+d^NN_kl(r), where d^NN_kl(r) is represented by an FNN (in this work always fixed to 16 hidden neurons) that inputs the system configuration r and outputs the symmetric matrix correction d^NN_kl. Notice that Ψ_PN is trivially a superset of Ψ_SN. 
    
    Although various architectures can be used, we adopt a three-layer fully-connected FNN for each of the functions a^NN_ki,σ and d^NN_kl (see Fig. <ref>).   The input layer has 2N neurons with neuron i (neuron i+N) outputting 1 if there is spin up (spin down) on site i and -1 otherwise, where N is the total system size. The hidden layer contains mN hidden neurons for constant m with Rectifier Linear Units (ReLU) <cit.> activation functions.   The output layer then contains O(N^2)   neurons specifying the values of the respective functions.  Gradients are computed in the standard way using variational Monte Carlo (see Supplementary Material  <cit.> Sec. II) which requires evaluating the derivative of the wave-function with respect to the weights and bias in the neural network.  Derivatives for FNN are typically taken using back-propagation. Because the wave-function is a determinant of a matrix generated by the neural-network output, we evaluate this full derivative by envisioning this determinant as an additional final layer of the neural network and then performing back-propagation including this layer. This ensures the cost of computing all the derivatives is of the same order as the evaluation of the wave-function (see Supplementary Material  <cit.> Sec. III).   Optimization is performed by stepping each parameter in the direction of the gradient with a random magnitude, which helps us avoid shallow local minima <cit.>, or by the RMSPROP method <cit.>. 
    
    
    
    The computational complexity of the NNB implemented with a single layer of O(mN) hidden neurons scales as O(mN^4) per sweep (i.e. after N electrons move) for forward and backward propagation and O(N^4) per sweep for the evaluation of the mean-field determinant.  This is similar to the scaling of standard backflow. 
    
    Explicit construction of standard backflow.—In this section, we provide an explicit construction which represents the standard backflow transformation in the form of Eq. (<ref>) written as a NNB. 
    
    In Eq. (<ref>),  η_ij,σ = t D_i H_j θ_|i-j|,σ=t n_i,↑n_i,↓ h_j,↑h_j,↓θ_|i-j|,σ, where h_j,σ=1-n_j,σ, θ_1,σ and θ_2,σ are the only non-zero variational parameters. We first demonstrate that η_ij,σ can be presented by a two layer neural network with input layer as (σ_1,...,σ_N,σ_N+1,...,σ_2N) where σ_i=2n_i,↑-1 and σ_i+N=2n_i,↓-1. By construction, n_i,σ and h_j,σ take value of 0 or 1 so that n_i,↑n_i,↓ h_j,↑h_j,↓ is 1 if and only if n_i,↑=n_i,↓=h_j,↑=h_j,↓=1. Therefore, t θ_|i-j|,σ D_i H_j = t θ_|i-j|,σ n_i,↑n_i,↓ h_j,↑h_j,↓ is equivalent to ReLU[t θ_|i-j|,σ (n_i,↑+n_i,↓+h_j,↑+h_j,↓-3)], which is the same as ReLU[t θ_|i-j|,σ (σ_i/2+σ_i+N/2-σ_j/2-σ_j+N/2-1)]. As a result, for each η_ij,σ, we associate it with a hidden neuron, such that the weights connecting it to σ_i,σ_i+N,σ_j,σ_j+N are tθ_|i-j|,σ/2,tθ_|i-j|,σ/2,-tθ_|i-j|,σ/2,-tθ_|i-j|,σ/2 respectively, the bias is -tθ_|i-j|,σ and the activation function is ReLU. In general, for more complicated backflow <cit.> with terms n_i,σh_i,-σ n_j,-σh_j,σ, n_i,σn_i,-σ n_j,-σh_j,σ and n_i,σh_i,-σ h_j,σh_j,-σ, where σ is the spin index, we can use more hidden neurons and represent it in the same way.
    
    After we have the neural network construction for the standard η_ij,σ, the term a^NN_ki =∑_jη_ij,σϕ_kσ(r_j,σ) in Eq. (<ref>) can be realized through an extra layer taking the outputs η_ij,σ to a neuron representing a_ki where the weight is given by the single particle orbital values ϕ_kσ(r_j,σ), there is no bias and the activation function is the identity.  This construction shows that the standard backflow parameterization is thus a subset of our three-layer NNB. 
    
    
    Results.—We have benchmarked the quality of our NNB on a number of systems including Hubbard models at various sizes and doping (all at U/t=8) as well as a frustrated magnet, the Heisenberg model on the Kagome lattice.  
    In the main text we focus primarily on the Hubbard model at n=0.875 (primarily on the 4 × 4 lattice) leaving the additional benchmarks as Supplementary Material  <cit.> Sec. IV., <cit.>)   
    The Hubbard Hamiltonian is
    
        H = -t ∑_iσ (c^†_iσ c_i+1σ + h.c.) +  ∑_i U n_i↑ n_i↓
    
    where we use U/t=8.  We compare the results to an optimized unrestricted (i.e. different single particle-orbitals for spin-up and spin-down) Slater Determinant (Ψ_S0) as well as a backflow BDG wave function (Ψ_PB) which transforms single particle orbitals of each spin by Eq. (<ref>). The formulation and the variational parameters of each wave function ansatz are summarized in Table. <ref>.  The parameters which aren't optimized, such as the initial set of orbitals {ϕ_kσ(r_i,σ)}  are obtained for  Ψ_PB and  Ψ_SN  by optimizing a restricted Slater-Jastrow wave-function  while Ψ_PN uses orbitals taken from the free hopping Hamiltonian (in practice the nature of the neural net allows for a direct change to the orbitals by altering the bias' on the final layer).  
    
    
    The relative error of the energy of NNB is 1.4% (and 0.66% after variance extrapolation (see Supplementary Material  <cit.> Sec. I) which is significantly better then the standard wave-functions (see Fig. <ref>(left)).  We examine the effect of the number of hidden neurons mN (see Fig. <ref>(right)). We find that at small hidden neuron number, Ψ_PN is much better than Ψ_SN but this advantage eventually largely disappears at large neuron number suggesting that a backflow parameterized with a small neural networks can compensate for the missing pairing in a Slater-determinant.  Surprisingly in the regime we've probed both NNB have energies linear with respect to 1/m in spite of the fact that in the m→∞ limit, they both must become exact as the FNN could simply put the exact amplitude Ψ(R) on one element of the diagonal <cit.>, one on the rest of the diagonal and zero everywhere else. 
    
    
    
    In addition, we have studied the NNB at 4 × L for L={4,8,12,16} comparing against the DMRG energy (PBC,open) for the 4 ×∞ system of -0.7659 ± 4 × 10^-5 per site <cit.>. We find our result very comparable (see Fig. <ref> (bottom)) to the DMRG result, especially for system sizes which are commensurate with wave-length 8 stripes<cit.>.
    
    We also investigate how the neural network backflow wave function affects the observables of our system. We see the expectation of double occupancy decreases as 1/m (see Supplementary Material  <cit.> Sec. I) and the spin and charge densities become significantly more symmetric as the number of hidden neurons increases (see  Fig.  <ref>).
    
    To understand the role of neural network in backflow transformation, we investigate how the neural network backflow transformation modify the orbitals. In Fig. <ref> (right), we notice although the backflow transformation on spin up orbitals and spin down orbitals are performed by two different neural networks, they produce similar backflow transformed orbitals and roughly preserves the symmetry of spin up and spin down for a given configuration. This is different from the optimized  unrestricted Slater Determinant Ψ_S0, which breaks the spin up and spin down symmetry significantly (see Fig. <ref> (Left)).
    
    One feature of using a NNB is the ability to alter the sign-structure of the wave-function.  Here we consider the amount the sign changes between Ψ_SN with 16 hidden neurons and Ψ_S0 by evaluating the integral 
    
        ∫ |Ψ_S0(x)|^2 sgn(Ψ_SN(x))  sgn(Ψ_S0(x)) dx/∫ |Ψ_S0(x)|^2 dx
    
    which is approximately 0.815 giving a 9% difference between the signs. 
    
    
    
    Furthermore, we open up the Ψ_SN neural network for m=8 and analyze the weight between the input layer and the hidden layer, which represents the features that the neural network learns from input. In Fig. <ref>, we plot these weights for both the spin-up and spin-down neural networks.  Interestingly the spin-up neural network primarily has large weights connected to the spin-down configurations while the spin-down neural network primarily has large weights connected to the spin-up configurations. This allows the neural network to introduce correlation between spin-up and spin-down configurations. Another observation is that more neurons tend to take large weight in negative bias, and small weight in positive bias.  
    
    
    
    
    Conclusion.—In this paper, we utilize the generality of artificial neural networks and the physical insight from backflow to develop a new class of wave function ansatz, the neural network backflow wave function, for strongly correlated Fermion systems on lattice. It achieves good performance for Hubbard model at nontrivial filling. We also show improvement on a kagome Heisenberg model in the supplement. While this work has focused on Fermion system on the lattice, the NNB is straightforward to generalize to frustrated spin systems as well as the continuum.  In the latter case, the input could be represented as a lexicographically ordered set of particle locations.  Our work provides a new approach toward combining machine learning methodology with dressed mean-field variational wave-functions which allows us to take simultaneous advantage of their respective strengths.  
    
        Acknowledgement.—This project is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the State of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This material is based upon work supported by the U.S. Department of Energy, Office of Science under Award Number FG02-12ER46875. Di acknowledges useful discussion with Ryan Levy, Dmitrii Kochkov, Eli Chertkov, Yusuke Nomura and Giuseppe Carleo. 
    
    
    
    § SUPPLEMENTARY MATERIALS
    
    
    § NUMERICAL RESULTS ON NEURAL NETWORK BACKFLOW OF STANDARD FORM
    
    
    We have shown that the standard backflow can be represented by a three layer neural network. In practice, we have also investigated the neural network in the following form, which is a direct representation of the standard backflow <cit.>
        ϕ_kσ^b(r_i,σ;r) =  ϕ_kσ + ∑_jη^NN_ij,σϕ_kσ(r_j,σ)
    
    where η^NN_ij,σ is represented by a three layer neural network similar to the one used for a^NN_ij,σ.
    
    
    
    
    
    We implement this backflow transformation in the Slater Determinant mean field, which is named as Ψ_SN-b. Notice that Ψ_SN-b can also be represented by Ψ_SN since we can define a^NN_ij,σ=∑_jη^NN_ij,σϕ_kσ(r_j,σ). In practice, Ψ_SN and Ψ_SN-b have similar performance in computing energy and doublon density (see Fig. <ref>, Fig. <ref> and Fig. <ref>).
    
    
    
    
    
    
    
    
    
    § OPTIMIZATION SCHEME
    
    
    Given a wave function ψ and Hamiltonian H, the energy E=⟨ψ|H ψ⟩/⟨ψ | ψ⟩. For any variational parameter t, to compute the derivative ∂ E/∂ t in quantum Monte Carlo, we define the local energy E_L=⟨r|H ψ⟩/ψ(r), where r is a system configuration. Then the derivative can be computed through Monte Carlo sampling.
    
    
        ∂ E/∂ t = 2 (⟨ E_L (1/ψ∂ψ/∂ t)⟩ -E ⟨1/ψ∂ψ/∂ t⟩)
    
    where ⟨⟩ is averaged over the Monte Carlo samples.
    
    To update the parameter t in the k+1 -th iteration, we use the sign of the gradient with a random magnitude as follows.
    
    
        t_k+1= t_k - αγ|∂ E/ ∂ t_k|/∂ E/ ∂ t_k
    
    where α is a random number in the range (0,1), γ is the step size, t_k and t_k+1 are the values for variational parameter t at the k-th and k+1-th iteration.  The scheme that alpha is random is inspired by the method described in <cit.>.  This method can help to get the optimization out of shallow local minima and is less sensitive to noise then standard gradient descent approaches. As a practical matter, we have optimized using this approach, gradient descent and stochastic reconfiguration and find that this approach is the best tradeoff between speed and quality of optimization. In practice, we choose γ=0.001, weights and bias' of the neural network are initialized from uniform distribution in the range of (-0.005,0.005).
    
    Notice that to represent Slater Determinant with NNB, one can set the bias of the output layer in NNB to be the orbitals while the other parameters to be zero. This observation can reduce the difficulty of optimization by initializing NNB with the optimized Slater Determinant, where the bias of the output layer is the optimized orbitals and the other parameters are small random numbers around zero. This trick is helpful for optimization in large systems.
    
    
    
    § BACKWARD PROPAGATION OF NEURAL NETWORK BACKFLOW
    
    
    From the optimization scheme, it is clear that we need to compute 1/ψ∂ψ/∂ t. In this section, we give the details on how to compute this quantity by backward propagation.
    
    In neural network algorithm, forward propagation refers to the evaluation of the network while backward propagation is essentially performing chain rule to compute derivatives. For forward propagation, the values in the i+1-th layer are computed from the values in the i-th layer under matrix operation and activation function, so that the output layer's values can be obtained from the input layer's values. For backward propagation, intuitively one can understand it as backward propagation of 'errors', where the 'errors' in the i-th layer are computed from the 'errors' in the i+1-th layer. More accurately, the 'errors', such as ^iδ^l_p in Eq. (<ref>), are related to the derivatives of the parameters. It indicates that the derivatives  of parameters in the i-th layer level can be computed from the derivatives of in the i+1-th layer level, which is essentially the chain rule of computing derivatives. A detailed description of applying backward propagation in NNB is given as follows.
    
    Consider backward propagation in a neural network with N layers. Denote the value in each layer by a set of vectors a^0, a^1, ..., a^N-1. Denote the value before activation in each layer (counting from the 2nd layer) by a set of vectors z^1,z^2,...,z^N-1. Denote the activation function in each layer (counting from the 2nd layer) by a set of functions f_1,f_2,...,f_N-1. Denote the weights by a set of weights (counting from the 2nd layer)  by a set of matrix w^1,w^2,...,w^N-1 and the bias by a set of vectors by a set of vector b^1,b^2,...,b^N-1. Then a^0 is the input configuration, for 1 < l ≤ N-1 , forward propagation reads
    
    
        z^l   = w^l a^l-1 + b^l
        
                a^l   = f_l (z^l)
    
    
    To calculate derivative, we use backward propagation. For the i-th output in the final layer, denote the error of each layer (counting from the 2nd layer) with respect to this output by a set of vectors ^iδ^1,^iδ^2,...,^iδ^N-1. Set ^iδ^N-1= f'_N-1(z^N-1_i). For 1 ≤ l ≤ N-1, backward propagation reads
    
    
        ^iδ^l-1   = ((w^l)^T  ^iδ^l ) * f'_l-1(z^l-1) 
        ∂ a^N-1_i/∂ b^l_p   = ^iδ^l_p 
        ∂ a^N-1_i/∂ w^l_pq   = a^l-1_q ⊗^iδ^l_p
    
    where * is the element-wise multiplication and the last equation uses outer product ⊗. 
    
    
    To transfer backward propagation to 1/ψ∂ψ/∂ t, we should view it as a special mapping on top of a^N-1. Therefore, we only need to choose a different ^iδ^N-1 and the rest will be the same as above. First, we consider the Slater-Determinant type neural network backflow.
    
    
        1/ψ_SD∂ψ_SD(r)/∂ t    =   1/det[M^SD,↑]∂[M^SD,↑]/∂ t
           +1/det[M^SD,↓]∂[M^SD,↓]/∂ t;
    
    where M_ik^SD,σ =ϕ^b_kσ(r_iσ), t is any variational parameter {w^l_pq,b^l_p}. Under the Einstein notation, the chain rule gives rise to, 
    
    
        1/det[M^SD,σ]∂[M^SD,σ]/∂ t
            =1/det[M^SD,σ]∂[M^SD,σ]/∂ϕ^b_kσ(r_iσ)∂ϕ^b_kσ(r_iσ)/∂ t
            =tr( (M^SD,σ)^-1∂[M^SD,σ]/∂ϕ^b_kσ(r_iσ) ) ∂ϕ^b_kσ(r_iσ)/∂ t
            = (M^SD,σ)^-1_ki∂ϕ^b_kσ(r_iσ)/∂ t
    
    
    For Ψ_SN, we have
    
        1/det[M^SD,σ]∂[M^SD,σ]/∂ t
            = (M^SD,σ)^-1_ki∂ϕ^b_kσ(r_iσ)/∂ a^NN_ki,σ(r)a^NN_ki,σ(r)/∂ t
            = (M^SD,σ)^-1_kia^NN_ki,σ(r)/∂ t
    
    Therefore, we set ^iδ^N-1_ki=(M^SD,σ)^-1_ki f'_N-1(z^N-1_ki) and perform the backward propagation.Notice that we use a matrix index ki for simplicity here. In practice, a matrix will be reshaped into an array for computation.
    
    
    For Ψ_SN-b, we have
    
        1/det[M^SD,σ]∂[M^SD,σ]/∂ t
            = (M^SD,σ)^-1_ki∂ϕ^b_kσ(r_iσ)/∂η^NN_ij,σ(r)η^NN_ij,σ(r)/∂ t
    
    We then set ^iδ^N-1_ij=(M^SD,σ)^-1_ki∂ϕ^b_kσ(r_iσ)/∂η^NN_ij,σ(r)  f'_N-1(z^N-1_ij) and perform the backward propagation. Notice that we use a matrix index ij for simplicity here. In practice, a matrix will be reshaped into an array for computation.
    
    
    Next, we consider Bogoliubov de Gennes type neural network backflow wave-functions Ψ_PN, 
    
    
        1/ψ_BDG∂ψ_BDG(r)/∂ t  =   1/det[Φ]∂[Φ]/∂ t
    
    where Φ_ij = ∑^N_k,l=1ϕ^b_k↑(r_i,↑) S_klϕ^b_l↓(r_j,↓), t is any variational parameter {w^l_pq,b^l_p} of a^NN_ij,σ and d^NN_kl.
    
    For the case that t is a parameter of a^NN_ki,↑, we rewrite Φ_ij = ∑^N_k,l=1ϕ^b_k↑(r_i,↑) R_kj, where R_kj= ∑^N_l S^NN_klϕ^b_l↓(r_j,↓), ϕ^b_k↑(r_i,↑) is given by Eq. (5) in the main body of the paper and S^NN_kl=S_kl+d^NN_kl. The chain rule gives rise to 
    
    
    
        1/det[Φ]∂[Φ]/∂ t
           =1/det[Φ]∂[Φ]/∂ϕ^b_k↑(r_i,↑)∂ϕ^b_k↑(r_i,↑)/∂ a_kl∂ a_kl/∂ t
           = 1/det[Φ]∂ det[Φ]/∂ϕ^b_k↑(r_i,↑)∂ a_kl/∂ t
           = Φ^-1_ji R_kj∂ a_kl/∂ t
    
    We then set ^iδ^N-1_ki=Φ^-1_ji R_kj  f'_N-1(z^N-1_ki) and perform the backward propagation.
    For t is a parameter of a^NN_ki,↓, similarly we define Q_il= ∑^N_k (ϕ^b_k↓)^T(r_i,↓) S^NN_kl and then set ^iδ^N-1_ki=Φ^-1_ij Q_jk  f'_N-1(z^N-1_ki).
    
    For t is a parameter of d^NN_kl, 
    
    
    
        1/det[Φ]∂[Φ]/∂ t
           =1/det[Φ]∂[Φ]/∂ S^NN_kl∂ S^NN_kl/∂ d_kl∂ d_kl/∂ t
           = 1/det[Φ]∂[Φ]/∂ S^NN_kl∂ d_kl/∂ t
           =  ϕ^b_k↑(r_i,↑) (Φ^-1)^T_ijϕ^b_l↓(r_j,↓)  ∂ d_kl/∂ t
    
    We set ^iδ^N-1_ki=ϕ^b_k↑(r_i,↑) (Φ^-1)^T_ijϕ^b_l↓(r_j,↓) f'_N-1(z^N-1_ki) and perform the backward propagation.
    
    
    
    
    § NEURAL NETWORK BACKFLOW ON VARIOUS SYSTEMS
    
    
    To further test the ability of NNB, we have optimized Ψ_SN in various systems: (i) 4 × 4 Hubbard model at U/t=8 with n=0.75 and n=1.0 (See Fig. <ref>, <ref>, <ref>). (ii) 8 × 4, 12 × 4, 16 × 4 Hubbard model with n=0.875 and U/t=8 (See Fig. <ref>,  <ref>, <ref>). (iii) 12 × 8 Hubbard model at U/t=8 with n=0.875 (See Fig. <ref>). (iv) 4 × 4 × 3 Kagome Heisenberg model with J=1 ((See Fig. <ref>)). The optimization has been implemented with the RMSPROP <cit.> method.
    
    For the 4 × 4 Hubbard model, we find that the smallest doping (n=0.75) does best with an error (variance extrapolated error) of 0.631% (0.233%) with the largest doping (n=1.0) doing the worst giving a 2.714% (1.745%) error (See Fig. <ref>, <ref>, <ref> for details). The results are compared with the exact diagonalization data <cit.>.
    
    We consider both the 16× 4 and 12 × 8 lattice at U=8 and n=0.875.  Because there are no exact answers at these system sizes, we compare against approximate (not necessarily variationally upper-bounded) AFQMC.  For our 16 × 4 system, using a number of hidden neurons  n_h={8,16,32,64}, we find the energy decreases from  -46.211 (for the optimized Slater determinant) to -47.745 (n_h=64). The variance extrapolated energy is approximately -49 which is reasonably close to the twist-averaged boundary condition AFQMC result of -49.088<cit.> (See Fig. <ref>). In the 12× 8 system we use n_h={8,16,32} with relative errors against the same system in AFQMC <cit.> of 6.3% for the optimized Slater determinant, 3.94% for the n_h=32 NNB and 0.655% after variance extrapolation  (See Fig. <ref>).  NNB significantly improves on the standard Slater-Jastrow methodology and the variance extrapolation result is competitive with  other techniques.
    
    
    Finding wave-functions for the ground state of the Heisenberg model on the Kagome lattice is challenging. The best per site energy of a  projected BDG result on a 4× 4× 3 lattice is approximately -0.4305 <cit.> compared to  the exact answer of -0.4387 <cit.>. We have implemented the neural network backflow on top of the Slater Determinant with  n_h = {8,16,32,64,128,256}. The associated energies ranges from -0.4311 (n_h=8) to -0.4339 (n_h=256), giving an improvement of the relative error of approximately 50%  (See Fig. <ref>).
    
    The results are summarized in Table <ref> and <ref>.
