{'fraction_non_alphanumeric': 0.03203479275504521, 'fraction_numerical': 0.009095971783924262, 'mean_word_length': 4.359579306424105, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Neural Network Quantum States with the Linear Method

Due to the strong correlations present in quantum systems, classical machine learning algorithms like stochastic gradient descent are often insufficient for the training of neural network quantum states (NQSs). These difficulties can be overcome by using physically inspired learning algorithm, the most prominent of which is the stochastic reconfiguration (SR) which mimics imaginary time evolution. Here we explore an alternative algorithms for the optimization of complex valued NQSs based on the linear method (LM), and present the explicit formulation in terms of complex valued parameters. Beyond the theoretical formulation, we present numerical evidence that the LM can be used successfully for the optimization of complex valued NQSs, to our knowledge for the first time. We compare the LM to the state-of-the-art SR algorithm and find that the LM requires up to an order of magnitude fewer iterations for convergence, albeit at a higher cost per epoch. We further demonstrate that the LM becomes the more efficient training algorithm whenever the cost of sampling is high. This advantage, however, comes at the price of a larger variance.

Due to the strong correlations present in quantum systems, classical machine learning algorithms like stochastic gradient descent are often insufficient for the training of neural network quantum states (NQSs). These difficulties can be overcome by using physically inspired learning algorithm, the most prominent of which is the stochastic reconfiguration (SR) which mimics imaginary time evolution. Here we explore an alternative algorithms for the optimization of complex valued NQSs based on the linear method (LM), and present the explicit formulation in terms of complex valued parameters. Beyond the theoretical formulation, we present numerical evidence that the LM can be used successfully for the optimization of complex valued NQSs, to our knowledge for the first time. We compare the LM to the state-of-the-art SR algorithm and find that the LM requires up to an order of magnitude fewer iterations for convergence, albeit at a higher cost per epoch. We further demonstrate that the LM becomes the more efficient training algorithm whenever the cost of sampling is high. This advantage, however, comes at the price of a larger variance.

I. INTRODUCTION
Neural network quantum states (NQSs) [1] constitute one of the most promising approaches for the description of quantum many-body systems. Since their introduction in 2017, they have been successfully applied to various lattice systems in one and two dimensions [2][3][4][5]. More recently they have also been applied to fermionic systems in second [6] and first quantization [7].
Although there is a diverse array of NQS architectures available, essentially all of the high precision ground state projection algorithms employ the Stochastic Reconfiguration (SR) [14] learning method. During the optimization, the SR updates the Ansatz wave function to mimic imaginary time evolution within the variational subspace, hence projecting onto the ground state. This corresponds to a stepwise rescaling of the gradient according to the quantum fisher matrix [15]. Thus, the learning is closely related to the natural gradient descent method [16].
There are two essential aspects in the evaluation of neural network architectures: (1) expressivity of the Ansatz and (2) efficiency of the learning scheme [17]. While a lot of effort has been put into the characterization of expressivity for different NQS architectures [18][19][20], the efficiency question has not received nearly as much attention. At present, systematic improvements to the training time can only be achieved by adjusting SR hyperparameters or by reducing the number of model * thorbenjan.frank@gmail.com parameters. The latter, however, results in a decrease of the expressivity of the NQS. The possibility of changing the training algorithm itself is currently inaccessible since at present, SR is the only reliable learning algorithm for NQSs. We take this as a motivation to implement a tailored learning algorithm for the optimization of NQSs, which is based on the Linear Method (LM), first developed in the context of Variational Monte Carlo (VMC) [21,22]. While, VMC typically uses real variational pa-rameters our results are, to the best of our knowledge, the first investigating the LM for wave functions with complex valued parameters.
In this paper we investigate the scaling of the computational cost when training NQSs with the LM and compare it with the SR. Such scaling analysis is ambiguous in the VMC setting which usually considers molecular systems, where it is difficult to distinguish effects arising form increasing the size of the simulation from those originating in the properties of the specific molecular orbital basis choice.
We find that the LM typically requires up to an order of magnitude fewer training epochs for learning the ground state, albeit at a higher cost per epoch. Thus, one faces a typical trade-off problem when searching for the more efficient learning algorithm. Indeed we find this problem to be non-trivial since it depends on the Hamiltonian under investigation as well as on the point in the Hamiltonian parameter space.
Our scaling analysis allows for the identification of a crossover point in which the choice for the more efficient training algorithm changes. Starting from there it can be shown that the location of this crossover point strongly depends on the computational cost of sampling. In this context sampling corresponds to the generation of samples {x k } w. r. t. the absolute square of the parametrized wave function |ψ θ (x)| 2 using Markov chain Monte Carlo (MCMC). More specifically we find that as long as sampling is easy (computationally cheap) the SR is more efficient. When sampling is expensive, however, the LM shows clear advantages with respect to computational training cost. These qualitative statements are substantiated in the main body of text.
Thus, we expect the LM to be a valuable tool for the training of NQSs whenever sampling requires large amounts of computational resources. We will illustrate the strength of the LM on the example of fermionic NQSs for simple chemistry Hamiltonians. There it has been shown that a tremendous amount of samples is necessary to generate sufficient statistics [6], which is due to the sharply peaked nature of the wave function around the Hartree-Fock state. Beyond that, our results are numerical evidence that the LM can be used for the learning of complex valued NQSs expanding the set of reliable and accurate learning algorithms.

A. Neural Network Quantum States
The basic idea of NQSs relies on parametrizing some quantum state vector in a basis {x} as such that the complex amplitude of the wave function ψ θ (x) is defined by the output of some neural network.
The parameter configuration θ * describing the ground state can be written as a minimization problem where E θ * is the global minimum of the variational energy H θ ≡ ψ θ |H|ψ θ in the n var -dimensional space.
As the exact evaluation of expectation values, and thus of the variational energy, scales exponentially in the system size, they are estimated by MCMC sampling is the classical expectation value for the current wave function ψ θ (x). Thus, the evaluation of E · · · M requires the generation of samples M = {x} distributed w. r. t. the absolute square of the current parametrized wave function |ψ θ (x)| 2 .
The value of the variational energy for some given parameter configuration can be calculated as where H loc (x) is the so called local energy Note that the local energy in principle depends on the point in the parameter space but its dependence on θ has been left out for notational convinience. The global minimum of equation (2) is found by iteratively updating the variational parameters θ where the update step ∆θ is calculated based on local approximations made to the variational wave function ψ θ and energy H θ . In the following we will review the two main learning algorithms that will be analyzed throughout this paper, the Stochastic Reconfiguration (SR) and the Linear Method (LM). Given that the LM is usually applied to real valued wave functions one has to carefully consider the complex case in order to make the LM applicable to NQSs.

B. Stochastic Reconfiguration
Originally introduced by Sorella in 2001 [14] for the optimization of Jastrow wave functions, the SR is currently the state-of-the-art method for the optimization of NQSs for quantum spin systems. The idea is to update the parameters in the variational wave function, such that it mimics the imaginary time evolution |ψθ ≈ e − H |ψ θ [18]. Beside its physical foundation, it is also closely related (and for positive wave functions equivalent) to the natural gradient descent [16], a learning algorithm originated in the machine learning community.
Within the SR each update step is calculated as where η is the so-called learning rate, S the quantum Fisher matrix and f the force vector For the readers convenience, a derivation of Eqns. (8,9) and of the the stochastic estimators is provided in appendices A and C. The quantity D k in the equations above describes the change in the wave function with the variational parameters such that ∂ θ k |ψ θ = D k |ψ θ , where ∂ θ k is the partial derivative w. r. t. to the k-th variational parameter θ k . Its corresponding stochastic estimator is the so-called log-derivative which is defined as The S matrix can be non-invertible due to vanishing eigenvalues. In order to ensure invertibility of S, a positive regularization constant a diag is added to its diagonal as where 1 is the identity.

C. Linear Method
The LM was originally introduced by Nightingale and Melik-Alaverdian [23] in 2001 and was later extended by Umrigar et al [22] to allow for the optimization of nonlinear parameters. At its core, the LM relies on the linear expansion of the explicitly normalized wave function, which lies in the self-plus tangent space Ω + [24], spanned by the current normalized wave function |ψ 0,θ and its tangent space Ω := |ψ | ψ|ψ 0,θ = 0 . Solving the projected Schrödinger equation in Ω + results in update steps which often outperform second order methods like the Newton-Raphson method, as they obey a strong zero variance property.
Lets start by considering the first order expansion of the explicitly normalized wave function |ψ 0,θ which can be written as where Φ is a global phase factor and |ψ 0,θ is the current normalized wave function with ∆θ 0 = 1. For k > 0 the derivative w. r. t. the k-th variational parameter is given as where D k are the log-derivatives (10). It can then be easily verified that the |ψ k,θ lie in the tangent space Ω of |ψ 0,θ for k > 0 such that |ψ θ+∆θ lies in the self-plus tangent space Ω + . While the orthogonality between |ψ 0,θ and |ψ k,θ is straightforward to show in the case of real parameters, the case of complex parameters generally results in a non-zero overlap between the current wave function and its derivatives, as pointed out by Motta et al [25]. However, up to first order, orthogonality can be restored for the complex case using phase optimization (see appendix B). With the linear expansion at hand, the idea of the LM relies on minimizing the linear energy approximation = min where we define H ≡ ψ θ+∆θ |H|ψ θ+∆θ and S ≡ ψ θ+∆θ |ψ θ+∆θ which are the so-called Hamilton and overlap matrix, respectively. The overlap matrix is fully expressed in terms of the log-derivatives D k where one identifies the quantum Fisher matrix (8) in the last expression. The entries of the Hamilton matrix also include the log-derivatives as well as additionally the Hamiltonian operator The expectation values can be estimated stochastically as in Eqn. (3) which is shown in detail for equations (16) - (22) in appendix C. The energy minimization as expressed by equation (15) is a quadratic problem in the parameter change ∆θ, where it will always hold that ψ θ+∆θ |ψ θ+∆θ ≥ 1 due to |ψ k,θ ∈ Ω for k > 0. Thus, equation (15) expresses the trade-off between energy minimization (making ψ θ+∆θ |H|ψ θ+∆θ as small as possible) and normalization conservation (keeping ψ θ+∆θ |ψ θ+∆θ as close to one as possible). The solution to (15) can be found as the vector v 0 ∈ C |θ| associated with the smallest real eigenvalue λ 0 of the generalized eigenvalue equation where c is a complex scalar. In order to ensure invertibility of S in a similar fashion as in equation (11) a shift is added to the diagonal of S apart from the first entry. The update within the LM is then given as

Stabilization by Regularization
Let |ψ be the new wave function that has been updated according to equation (24). Then its normalized form will be a good approximation of the subspace eigenfunction |ψ θ+∆θ ∈ Ω + , as long as the update step taken in the variational space is small enough [24]. Sufficiently small update steps can be ensured by adding a positive, real constant κ to the diagonal of the H matrix (except for the first entry) [26]. This regularization scheme plays the same role as the trust radius in the Newton-Raphson optimization and can be interpreted as damping of the update step, which means increasing κ reduces the length of δλ. The optimal choice of κ is by far not obvious and depends on the current point in the optimization as well as the system under investigation as it determines the energy landscape. A common way is to solve the eigenvalue problem defined in equation (23) for three different values of κ with κ n = κ 0 × 10 n . Afterwards, one chooses the update step that yields the lowest variational energy, which can be estimated over the correlated sampling technique described in [27].

D. Algorithmic Outline
For the SR and the LM, each training iteration consists of the following steps: (1) Generate n s samples with respect to |ψ θ (x)| 2 using the MCMC sampling.
(2) Calculate the stochastic estimates E · · · M for the quantum expectation values · · · that appear in the update matrices. Solve the linear set of equations for ∆θ.
Steps (1) -(3) are repeated until the energy E θ is considered optimal. A detailed analysis of the computational complexity can be found in the appendix D.

E. RBM Quantum States
Throughout, we will consider RBM quantum states, parametrized as |ψ θ = x ψ θ (x) |x where the complex amplitude associated with the basis states is given by the exponential family The binary vectors x and y of length N and M represent the visible and the hidden state, respectively. In the context of RBM quantum states the visible state is taken to be the basis state x. Each connection between hidden y i and visible node x j has an associated weight w ij and the vectors a ∈ C N and b ∈ C M describe the biases attached to visible and hidden nodes, respectively. We denote the whole set of variational parameters in the RBM quantum state as θ = {a, b, vec(w)} with weight matrix w ∈ C N ×M . The ratio between number of hidden and visible units α = M/N is commonly referred to as hidden unit density and enters the total number of variational parameters via n var = αN 2 + αN scaling quadratically in the system size N . For systems exhibiting full translational symmetrythe effective number of variational parameters can be reduced by one order of magnitude to n var ∝ O(αN ) (for more details see [1]). The optimization of RBM quantum states with the SR and the LM is computationally efficient when the following requirements hold [15]: (1) The computation of the operators D k (x), H loc (x) (and H loc,k (x) for the LM) is efficient for every state x.
(2) Efficient sampling with respect to ψ θ (x) is possible for any values of θ. This means the updates in a Monte Carlo chain can be computed efficiently, as the ratio of wave functions ψ θ (x) and ψ θ (x ) can be computed efficiently for any x and x .
(3) The Markov chain converges to the desired state in sub-polynomial time.
Although requirement (3) can hardly ever be guaranteed one observes in practice that it is often satisfied, as long as the Hamiltonian is free of frustration (note however recent proposals involving auto-regressive learning strategies, that guarantee efficient sampling [12]). Indeed it should be noted, that the requirements above are sufficient for the learning algorithms presented in this paper, but might be altered if one is using other, e.g. second order methods, that additionally rely on the efficient calculation of second order derivatives D k,k (x).

F. Spin Models
Throughout this paper, the LM and SR will be compared for the transverse field Ising (TFI) model as well as for the J1J2 model. In this context X, Y, Z will denote the Pauli operators with corresponding eigenvalues x, y, z = ±1 and periodic boundary conditions are assumed for all models. Beyond that, we apply local basis transformations to make the TFI model stoquastic and the J1J2 model less frustrated. This has been shown to improve the performance when training NQSs [28].

Transverse Field Ising Model
The Hamiltonian of the TFI model is given as where · , · denotes nearest neighbours and h the strength of an external field, which is applied perpendicular to the z axis. At h = 1 the TFI model has a quantum critical point while it reduces to the classical Ising model for h = 0. In 1D, the model has an exact solution by mapping it to free Fermions [29].

J1J2 Model
The second spin model that will be investigated is the J1J2 model which includes next nearest neighbour interactions. Its Hamiltonian is parametrized by the nearest and the next nearest neighbour couplings J 1 and J 2 : where · , · indicates next nearest neighbours. Throughout the paper, we fix J 2 = 1 and restrict J 2 ≥ 0 resulting in anti-ferromagnetic interactions such that the total magnetization of the ground state is zero, allowing for spin exchange sampling. For J 2 = 0 one recovers the anti-ferromagnetic Heisenberg model and the system is frustration free with long-range Neel order [30,31]. In 1D the model is exactly solvable at the Heisenberg point (J 2 = 0) as well as for J 2 = 0.5, where it becomes the Majumdar-Gosh (MG) model [32] and undergoes a phase transition to a frustrated phase. A second phase transition from a gapless spin-fluid regime to a gapped dimer regime is found within the ferromagnetic regime at J 2 ≈ 0.24 [33,34]. The model is known to be non-stoquastic for J 2 = 0, and poses significant challenges for VMC [28].

Chemistry Hamiltonians
Finally, we consider molecular Hamiltonians in secondquantized form where a † i and a j are fermionic creation and annihilation operators and Q ij and R ijkl are one and two-body integrals. Interacting fermionic Hamiltonians can be mapped to spin systems by using e.g. the Jordan-Wigner (JW) or the Bravyi-Kitaev (BK) transformation [35,36]. In the context of chemistry, NQSs using the JW transformation have been reported to yield more accurate VMC energies than those using the BK transformation [6]. Thus, we will concentrate on the JW transformation, which is formally given as where 2P ± j = X j ± iY j . Applying the transformation allows for mapping a Hamiltonian of the form given in equation (29) to an interacting spin Hamiltonian such that we can train RBM quantum states to represent their ground state. A. Scaling Analysis

Training Epochs
The learning curves for the TFI model in figure 1 suggest that training with the LM can require up to an order of magnitude fewer epochs than the SR to achieve comparable accuracies. This is illustrated more extensively in figure 3 where the number of epochs for convergence is plotted for various values of h and system sizes N .
We see that the LM convergences in a number of update steps which is typically an order of magnitude smaller than for SR. This is true of other second order methods, like the Newton-Raphson method. However, the LM is generally cheaper as second order derivatives do not need to be computed. Moreover one observes a qualitatively different behavior of n conv as a function of the external field h for SR and LM when approaching the quantum critical point h = 1. The number of epochs for convergence reaches a maximum near the critical point for SR, whereas the LM seems to be quite insensitive. This might indicate that in addition to faster convergence properties the LM is less sensitive w. r. t. quantum criticality.

Update Algorithm
Looking at the bare number of training epochs, however, is only half the truth in neural network training. The final quantity of interest is the clock time required to find a proper approximation for the ground state wave function. Indeed the LM has a much larger per epoch cost, which can in certain circumstances cancel the convergence speedup. To approach this issue we define the time for a single epoch as the sum of the time required for sampling (algorithmic step (1)) and the time required for the calculation of the update step (algorithmic step (2) and (3))  Each tile is colored according to the difference in the total update time ∆Tu = T sr u −T lm u . The total update time describes the clock time until convergence, which is required for the calculation of the update step within the SR and the LM, respectively (see eq. (33)). Essentially, this difference describes the training time difference between SR and LM when sampling is not taken into account. Blue tiles (∆Tu > 0) correspond to regions where the total update time for the LM until convergence is smaller and vice versa. (b) Runtime Phase Diagram for the Total Training Time in the TFI model: Here the same situation as in (a) with the time for sampling taken into account is shown. The difference in total training time is defined as ∆T = T sr − T lm , where the total training time is the clock time required for the training of the NQS (see eq. 34). Thus, it allows to identify the more efficient training algorithm. Again blue tiles are regions in which the LM is more efficient.
T u . It is defined as where t (i) u is the time for the update calculation at training epoch i. Thus the total update time T u scales linearly in the number of epochs. Figure 4.a shows what we call the runtime phase diagram for the total update time. It serves the purpose of globally identifying regions in which either the LM or the SR requires less total update time T u . Therefore, each tile (corresponding to a combination of external field h and system size N ) is colored according to the difference ∆T u = T sr u − T lm u such that ∆T u > 0 (blue) corresponds to regions in which the LM is faster w. r. t. the total update time. We find that the total update time for SR is shorter for most system sizes N and external field values h where one observes the difference to increase with increasing system size indicating asymptotic advantage. Only for an external field value of h = 0.9 do we see superior LM performance. Note that that each row in the runtime phase diagram corresponds to a plot like the one in figure 3 where the number of training epochs is replaced by the total update time. Thus, when looking only at the time that is required for the update algorithm itself, the SR is found to be widely superior in the asymptotic limit.

Total Training
So far, the time for sampling has been disregarded in the scaling analysis in order to make a clear statement about the scaling relation of LM and SR as stand alone components. Clearly, the quantity of ultimate interest, however, is the total time that is required for the training of NQSs Our simulations show clear evidence that, when taking sampling into account, the more efficient algorithm shifts from the SR to the LM (compare figure 4.a and 4.b). For the TFI model the LM converges to the ground state in less training time than the state-of-the art SR algorithm across most of the N −h-space ( figure 4.b). Moreover, we find the training time difference to increase when going to larger system sizes suggesting an asymptotic training time advantage for the LM.

B. Efficiency Regimes
The scaling analysis from above suggests that the sampling cost is the determining factor separating the efficiency of the SR and LM methods. This is clearly seen As the ground state of the J1J2 model lies in the M 0 z symmetry sector, MCMC exchange sampling can be used, which is computationally more efficient than single spin flip sampling, used for the TFI model (green and blue triangles in figure 5). A vastly more expensive sampling cost arises in the context of molecular Hamiltonians (red triangles in figure 5), the reason for which is twofold: First, particle number and spin conservation has to be respected throughout the sampling. To accomplish that we use the MetropolisHamiltonian sampler from the package NetKet [37]. Second, the total number of samples required is of the order O(∆k × 10 4 ), where ∆k = 10N is the downsampling interval. This is due to the sharply peaked nature of the absolute square of the wave function around the Hartree-Fock state [6], which results in low acceptance probabilities during sampling. Based on the number of training epochs n conv and the update time readily calculated by solving the equality We interpret the gap between sampling time per epoch and the transition line t trans s as a measure for how much more efficient either of the training algorithms is for a model. Within our numerical simulations we found the SR to be the more efficient training algorithm for the sampling cheap J1J2 model, as the sampling time per epoch t s lies below the transition line (green triangles and green solid line in figure 5). When sampling becomes more expensive as is the case for the TFI model, we find the LM to be more training efficient (blue triangles above blue solid line). As the gap between t s and t trans s grows in the number of sites we expect the LM to be asymptotically more training efficient than the SR. This is indeed consistent with our findings from the runtime phase diagram in the section above (see also figure 4.b).
For molecular NQSs, our numerical results indicate a strong training time advantage for the LM, where N = 4 and N = 12 correspond to the qubit representations (see eq. (30) and (31)) of the H 2 and the LiH molecule, respectively. Having the transition line close to zero furthermore indicates that the LM and the SR are nearly equally efficient when sampling is not taken into account. This is due to the fact, that the number of training epochs for molecular Hamiltonians with the LM is more than an order of magnitude smaller than with the SR.
Finally, it is important to note that any improvements on the sampling strategies for specific models can tip the balance in favor of the SR method.

A. Learning Reliability
In order to determine the learning reliability, we consider a VMC run to be converged if meaning if the relative error is below some threshold b.
If not stated otherwise we choose b = = 2 × 10 −3 . The relative error is given as where E 0 is the exact ground state energy. The convergence ratio is calculated as the fraction where N * vmc is the number of VMC runs that are labeled as converged according to eq. (37) and the total number of VMC runs is denoted as N vmc . We find high convergence reliability for h ≥ 0.3 and most values of J 2 , but also observe regions in the model parameter space where c r ≤ 0.25 ( figure 6). For the TFI chain c r is particularly low close to the classical point h = 0 and increases rapidly when moving away from it. Importantly the LM shows high reliability at and around the quantum critical point h = 1. For the J1J2 chain c r drops if one is deep in the gapped dimer regime at J 2 = 0.45. Despite similar convergence ratios for h ≤ 0.3 and J = 0.45, the reasons for the low convergence ratios of TFI and J1J2 model, however, originate from different causes. VMC runs that fail in the gapped dimmer regime usually come On the contrary, the non-convergent runs in the TFI get mostly stuck in local minima of the energy landscape with VMC energies far away from the real ground state energy. For that reason repeated convergence is only achieved when the threshold for the relative error is increased to b = 0.5 (dashed line with square markers in figure 6.a). More generally, within the J1J2 model we observe that slightly increasing the threshold also increases the value of c r to one for J 2 ≤ 0.3. This is consistent with our observation of larger fluctuations in the energy within the J1J2 chain than in the TFI chain (for h ≥ 0.5). Compared to the SR (red lines in figure 6) we find that the LM is less reliable for b = .
In particular, the SR shows also high reliability in regions in which the LM suffers from convergence ratios c r ≤ 0.25.

B. Accuracy
Our results in figure 7 show that the achieved accuracies with the LM strongly depend on the Hamiltonian parameters. In the TFI model we observe constant accuracies with a relative error rel ∼ 10 −4 . However, for h = 0.1 we find the LM to yield less accurate results rel > 10 −3 which is in line with our observations that the learning reliability is small for h close to the classical point ( figure 6). Thus, this regime can be considered to be particularly hard for the LM. For the J1J2 model, we find the results for the ground state energy to be mostly less accurate than the ones observed in the TFI chain by approximately one order of magnitude. For J 2 ≤ 0.3 the relative error is constant and increases when being deep in the gapped dimer regime. As for the TFI model this is consistent with a small convergence reliability c r at J 2 = 0.45. Compared to the SR, one finds the LM to be several orders of magnitude less accurate at h = 2. However, close to the quantum critical point h = 1 as well as in the gapped dimer regime of the J1J2 model, we find the accuracy for the SR and the LM to be of similar order.

V. DISCUSSION
In this paper, we explored the LM as a learning algorithm for complex NQSs. We presented numerical evidence that the LM can be used for the learning of complex valued NQSs, to our knowledge for the first time. At the same time, we found the LM to deliver accurate representations for the ground state of the TFI and the J1J2 model on a one dimensional chain as well as for molecular Hamiltonians in a spin basis. We observed the LM to converge in an order of magnitude fewer epochs than the state-of-the-art SR algorithm. As a consequence we found the LM to be the more time efficient training algorithm, whenever MCMC sampling becomes expensive in the NQS training. The particular strength of the LM is most evident on the example of molecular Hamiltonians, where we found it to be clearly advantageous compared to the SR w. r. t. training time. Thus, we expect the LM to be a valuable tool for the optimization of NQSs whenever sampling is expensive.
On the flip side, we observed a significantly larger vari-ance for optimization with the LM which resulted in less accurate and less reliable results compared to the SR. Moreover, the LM showed a fragile nature during training when the optimization lacked an appropriate sign structure to start with, such that we encoded it by using the Marshall sign rule. This is indeed consistent with prior observations made during the optimization of NQSs [3,28], but we found it to be more severe for the LM than for the SR. Consequently, a deeper investigation regarding better phase learning with the LM might be a topic for future research. With its origin in the chemistry community and thus usually applied to problems in first quantization, it might be an interesting direction to investigate the LM in the context of recently introduced NQSs in first quantization basis [7]. Also in the context of chemical systems, the LM has been used for the direct optimization of excited states [39], such that one could test an LM based optimization for excited NQSs similar to approaches with the SR [2].
where we assume that |ψ 0,θ is holomorphic. Using Wirtinger derivatives ∂ θ k and ∂ θ * k the first order expansion of | ψ 0,θ can be written as (see below for full details) where the change with respect to the k-th variational parameter is given as In equation (B3), D k = ψ 0,θ |D k | ψ 0,θ denotes the expectation value of the log-derivatives (see eq. 10 main text). Consequently, in the case of real parameters the current wave function and its derivatives build a self-plus tangent space as the orthogonality appears quite naturally as tangent space of a sphere in real space. Here we assumed that real network parameters also result in real variational wave functions, which is indeed the case for NQSs but e.g. does not hold for variational quantum eigensolvers (VQEs). On the contrary, complex parameters lead to a purely imaginary overlap between the current wave function and its derivatives indicating that | ψ k,θ / ∈ Ω as also suggested in [25]. Consequently, the first order expansion of the normalized wave function | ψ θ+∆θ does not lie in Ω + for the complex case in general. We will now show that it is possible to construct a first order expansion |ψ θ+∆θ , which lies in Ω + . Starting point is the linear expansion given in equation (B2) where we have introduced with |ψ 0,θ ≡ | ψ 0,θ . It can easily be verified that |ψ k,θ ∈ Ω, by multiplying ψ 0,θ | from the right to (B6) and realizing that the overlap is indeed zero. Adding a term quadratic in δλ k to equation (B5), such that the pre-factors become the first order Taylor expansion of the exponential function, allows to write [40] |ψ λ+∆θ k = 1 + i Im D k ∆θ k |ψ 0,θ + 1 + i Im D k ∆θ k ∆θ k |ψ k,θ + O |∆θ k | 2 = exp iφ k |ψ 0,θ + ∆θ k |ψ k,θ with φ k = Im D k ∆θ k . The general expression when several parameters are changed can easily be deduced from equation (B7) and is given by where Φ = k φ k , is a global phase factor. Choosing Φ such that the pre-factor becomes one (as it will not affect the following derivations), the first order expansion of the normalized wave function can be written as where |ψ n,θ ∈ Ω for n > 0 by construction and consequently |ψ θ+∆θ ∈ Ω + . It should be noted, that the construction of |ψ θ+∆θ ∈ Ω + corresponds to the heuristic re-scaling procedure outlined in [22] with ξ = 1, which allows for the update of non-linear parameters. This underlines the importance of the derivatives to lie in Ω, which thus is considered to be crucial for the stability of the LM [41].

Log-derivatives
The scaling of the log-derivatives itself scales as O n var (by definition from above). Consequently, the stochastic estimation of all log-derivatives scales as O n s × n var .

Local Energy Derivatives
Recalling that the local energy derivatives are given as where H k (x) is given as Consequently, the calculation of all local energy derivatives for given x scales as O N × n var × M . The stochastic estimate scales then as O n s × N × n var × M .

Stochastic Reconfiguration
Within the SR, one has to construct the quantum Fisher matrix S and the force vector f which have the dimensions (n var × n var ) and (n var × 1), respectively. Consequently, their construction scales as O n 2 var and O n var . The update step is found as solution to a linear set of equations build by S and f ( §A). Its solution can found by building the inverse S −1 , which scales as O n 3 var .

Iterative Solver
In order to speed things up, one can make use of an iterative algorithm proposed in [44], which erases the need for explicitly building the S matrix as well as calculating its inverse. Instead, the linear set of equations is solved by an iterative conjugate gradient method, which requires a computational cost of O N var per step, due to the product structure of S [1]. Usually, the number of CG iterations is much smaller than the number of variational parameters such that using an iterative CG scheme greatly reduces the computational cost of the SR algorithm.

Linear Method
The construction of the Hamilton H and overlap matrix S scales as O n 2 var thus is in line with the scaling for the construction of S in the SR. However, as the LM requires calculation of the local energy derivatives (see (D3)) the calculations of the matrix entries itself is computationally more expensive than in the SR. The update step is found as solution to the generalized eigenvalue problem of dimension n var × n var , which generally scales as O n 3 var (build the inverse on one side and solve the resulting regular eigenvalue equation with the Schur decomposition). The scaling for the eigenvalue problem can be reduced by using the Implicitly Restarted Arnoldi Method (IRAM).

Implicitly Restarted Arnoldi Method
IRAM is particularly strong when one is interested in the outer eigenvalues/vectors of a system, which is the case in the LM (eigenvalue/vector with smallest real part). For an extensive introduction to theoretical and numerical details the reader is referred to [45] and [46]. As it is an iterative solver the actual gain in computational cost is hard to estimate and strongly depends on the hyper parameters and the structure of the matrices. To our knowledge finding appropriate upper bounds for the run time of IRAM is still an open question.

Network Architecture
For simulating the TFI and the J1J2 model we used symmetry encoding RBM quantum states with hidden unit density of α = 2.
The chemistry Hamiltonian simulations were performed without the symmetry encoding RBM architecture and also with a hidden unit density of α = 2.

Training Hyperparameters
The following table lists the training hyperparameters used for the data shown in the different figures.