Is Variational Monte Carlo Robust? Sharp Moment Thresholds and Heavy-tailed Stochastic Optimization

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Is Variational Monte Carlo Robust? Sharp Moment Thresholds and Heavy-tailed Stochastic Optimization

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Is Variational Monte Carlo Robust? Sharp Moment Thresholds and Heavy-tailed Stochastic Optimization Philipp Grohs1,2 and Davide Nobile1

arXiv:2606.26009v1 [cs.LG] 24 Jun 2026

Faculty of Mathematics, University of Vienna [email protected] 2 RICAM, Austrian Academy of Sciences [email protected] June 25, 2026 Abstract

Variational Monte Carlo (VMC) is a central algorithm in electronic structure theory and has gained renewed importance through modern neural-network ansätze such as FermiNet. At its core, VMC seeks ground states by minimizing the Rayleigh quotient by stochastic optimization. In this work, we show that the resulting stochastic optimization problem is intrinsically governed by the nodal geometry of the underlying wave function. More precisely, we establish that properties of the nodal set determine the integrability of the local energy and gradient estimators that drive VMC. For broad and practically relevant ansatz classes, including Slater-Jastrow wave functions with variable-exponent Slater-type orbitals, we prove that these estimators are generically heavy-tailed and fail to admit higher moments. At the same time, for general analytic ansätze, we prove weak moment bounds for the relevant estimators and identify precise low-moment regimes, showing how generic and degenerate nodal structures lead to different integrability thresholds. Building on this analysis, we introduce a new robust variant of VMC – coined PS-Clip-VMC – which is based on clipping both the local energy and the gradient random variable. We prove that PS-Clip-VMC converges both in expectation and with high probability in the weak moment regime of VMC. Preliminary experiments for training FermiNet on Atoms with up to 18 electrons suggest that PS-Clip-VMC is significantly more robust than standard methods.

Introduction

1.1

Motivation

A central computational problem in quantum chemistry concerns the approximation of ground states and ground state energies of many-particle Hamiltonians. If H : H 2 (Ω) → L2 (Ω) is a self-adjoint Hamiltonian bounded from below, its ground-state energy E0 is characterized by the Rayleigh–Ritz variational principle E0 :=

inf ψ∈d(H)\{0}

⟨Hψ, ψ⟩ , ∥ψ∥2L2

(1)

with d(H) the form-domain of H; typically we have d(H) = H 1 (Ω) [42]. If the infimum is attained, then every minimizer is an eigenfunction associated with the lowest eigenvalue of H – a so-called ground state. In many-particle systems we usually have that Ω = Rd×N (with d the particle dimension and N the number of PN PN particles) and (in dimensionless units) H = −∆ + i,j=1 V (xi , xj ) + i=1 Vext (xi ), where V : Rd × Rd → R describes particle interactions and Vext : Rd → R describes an external potential. If the particles are indistinguishable and fermionic, there is an additional restriction, namely ψ ∈ Ha1 (Rd×N ), where we denote by Ha1 (Rd×N ) the space of all functions ψ ∈ H 1 (R3×N ) with ψ(xπ(1) , . . . , xπ(N ) ) = sgn(π)ψ(x1 , . . . , xN ) 1

for all permutations π of {1, . . . , N }, all (x1 , . . . , xN ) ∈ Rd×N and sgn(π) denoting the signature of the permutation π. We generally refer to multi-particle systems restricted to the space of antisymmetric functions as fermionic multi-particle systems and note that for such systems, the problem (1) amounts to computing E0 :=

inf

ψ∈Ha1 (R3×N )\{0}

⟨Hψ, ψ⟩ . ∥ψ∥2L2

(2)

A particularly important example is given by the (Born-Oppenheimer) multi-electron Hamiltonian HBO := −∆ +

(3)

I=1

which describes electronic properties of a molecule with M nuclei (R1 , . . . , RM ) ∈ R3×M and charges (Z1 , . . . , ZM ) ∈ NM . The multi-electron Hamiltonian acts on the electronic wavefunction ψ ∈ H 2 (R3×N ) which depends on the n electron coordinates (x1 , . . . , xN ) ∈ R3×N [18, 41]. Since electrons are indistinguishable fermionic particles, (and ignoring spin) the electronic wavefunction is antisymmetric with respect to permutations of the electron coordinates. Finding efficient algorithms to accurately compute the groundstate energies and ground states of (3) (along with excited states corresponding to higher eigenvalues and eigenvectors) is one of the key open problems in computational chemistry. Its solution would enable the simulation and prediction of all nonrelativistic properties of molecules from first principles without resorting to expensive and time-consuming experiments. In the case of the Born-Oppenheimer multi-electron Hamiltonian (3), the solution of (2) poses formidable challenges: First, the dimension of the problem scales linearly in the number of electrons. This means that even for small molecules with only tens or hundreds of electrons, solving (2) requires solving a PDE eigenvalue problem on a computational domain with hundreds or even thousands of dimensions. Second, many chemical applications require exceedingly accurate approximations of E0 to within a few mHa (millihartree) while the absolute energy E0 may be of the order of several hundred Hartree. In other words, practically useful computations often need to approximate E0 to within a relative accuracy of 6 − 7 digits [11]. In recent years, so-called variational monte-carlo (VMC) algorithms have emerged as a new state of the art in terms of accuracy, see for example [31, 19, 34, 37, 35, 36, 38, 12, 11, 13, 20, 8, 9, 10]. These algorithms build on a parametric function class Ψ = {ψθ : θ ∈ Θ} ⊂ Ha1 (R3×N ) \ {0} and Θ ⊂ RP (for example consisting of neural networks) and approximate the true ground state and ground-state energy by minimizing the loss ⟨Hψθ , ψθ ⟩ L(θ) := . (4) ∥ψθ ∥2L2 Since evaluating L(θ) or ∇θ L(θ) requires the calculation of high-dimensional integrals, stochastic representations and corresponding stochastic optimization algorithms have to be used: with Xθ ∼

Hψθ (Xθ ) ∇θ ψθ (Xθ ) |ψθ (x1 , . . . , xN )|2 , Eθ := , Wθ := and Zθ := Eθ · Wθ − E [Eθ ] · E [Wθ ] ∥ψθ ∥2L2 ψθ (Xθ ) ψθ (Xθ )

we have that L(θ) = E[Eθ ]

and

∇θ L(θ) = 2E[Zθ ].

VMC algorithms leverage these stochastic representations by first sampling from Xθ via MCMC and then estimating ∇L(θ) from suitable finite sample estimators of E[Zθ ] to compute parameter updates [43, 5]. The success of this approach hinges on the quality with which the true gradient can be approximated by its respective stochastic approximation. Indeed, if the random variables Eθ and Zθ satisfy sufficient moment assumptions, [1] recently proved convergence of VMC to stationary points of L.

1.2

Contributions

The main contributions of this paper follow three threads which we now summarize.

1.2.1

Non-Existence of Moments

Our first main result uncovers yet another challenge in VMC: the moment assumptions of [1] are in general not satisfied and both Eθ and Zθ are severely heavy-tailed. More precisely, in Theorem 2.4 we prove that even for the simple case of Ψ consisting of any (sufficiently expressive) Slater-Jastrow Ansatz with SlaterType orbitals [34] the set of parameters with heavy tails on Eθ and Zθ is residual (in fact, we show that it contains a dense open subset). For h1 , . . . , hN ∈ C(R3 ) and (x1 , . . . , xN ) ∈ R3×N we write   h1 (x1 ) . . . hN (x1 )   .. .. |h1 , . . . , hN ⟩(x1 , . . . , xN ) := det  ... . . . h1 (xN )

...

hN (xN )

Then, we have the following result. Theorem 1.1 (Colloquial Version of Theorem 2.4). With Yl,m denoting the spherical harmonics and α ∈ 2 [0, ∞)I , c = (cl,m,k,i )l=0,...,L, m=−l,...,l, k=0,...,K−1, i=1,...,I ∈ R(L+1) KI denote X x |x|l+k e−αi |x| , h(α,c) (x) := cl,m,k,i Yl,m |x| l,m,k,i

the corresponding Slater-Type orbital and let Jβ be a symmetric and positive Jastrow factor, parametrized by β ∈ RJ . 2 For θ := β, ((αba , cba ))a=1,...,N, b=1,...,B ∈ RJ × [0, ∞)N BI × RN B(L+1) KL =: Θ denote ! B X ψθ (x) := Jβ (x) · |h(αb1 ,cb1 ) , . . . , h(αbN ,cbN ) ⟩(x) . b=1

Then, whenever L ≥ 1 and KI ≥ 3N − 3, the set n θ ∈ Θ : E[|Eθ |3 ] = ∞

∧

o E[|Zθ |3/2 ] = ∞

is residual. Our result in particular shows that the gradient random variable Zθ generically does not possess finite variance, rendering standard stochastic approximations insufficient. The main difficulty in the proof of Theorem 2.4 is to show that every neighbourhood of every θ in Θ contains a parameter µ such that ψµ has a regular zero (x1 , . . . , xn ) ∈ R3×n with ∆ψµ (x1 , . . . , xn ) ̸= 0 and ∇θ ψµ (x1 , . . . , xn ) ̸= 0, see Lemma 2.11. Hψ Having established this, the lack of moments follows from asymptotically expanding the local energy ψµµ and gradient random variable

Hψµ ∇θ ψµ ψµ · ψµ

around this zero, see Lemma 2.5.

Remark 1.2. Theorem 2.4 holds for any Hamiltonian. In particular, the lack of moments is not caused by singular potentials such as (3). Moreover it is possible to extend the result to general Slater-Jastrow Ansätze whenever the orbitals are chosen from a sufficiently expressive space, see Lemma 2.11. Remark 1.3. At present, our main result applies to Slater-Jastrow Ansätze, but in recent years more general representations have become popular. For example, in [19, 34] the Slater-Jastrow Ansatz is generalized by applying ψθ to correlated variables Bθ (x) through a parametrized backflow transform. A more general ansatz is given by FermiNet [31] and subsequent work [39, 37, 35, 36, 38, 12, 11, 13, 20, 8, 9, 10]. We expect that the conclusion of Theorem 2.4 remains true in these more general models using similar methods as in the proof of Theorem 2.4. We leave this to future work. 1.2.2

Existence of Moments

Having established a generic lack of high moments, our second main contribution consists of sharp lower bounds on the number of moments of Zθ : For real-analytic parametric function classes Ψ, we show in Theorem 3.2 that (subject to additional compactness assumptions) a minimal number of moments always exists: 3

Theorem 1.4 (Colloquial Version of Theorem 3.2). Let K ⊂ Ω compact with smooth boundary, and let η ∈ C ∞ (Ω) be any smooth cutoff function with supp(η) = K. Then, for any parametrized family Θ ∋ θ 7→ ψθ := ηφθ with φθ real-analytic and Θ compact, it holds that sup{p ∈ (0, ∞) : sup E[|Zθ |p ] < ∞} ≥ θ∈Θ

5 . 4

The proof of Theorem 3.2 makes use of the Weierstrass preparation theorem [33, Section I.3.3], along with an auxiliary non-degeneracy result in Lemma 3.6 to carefully control the local structure of nodal sets of ψθ . A few remarks are in order. Remark 1.5. By considering the function x 7→ x21 it is easy to see that the threshold p = 45 is sharp. Remark 1.6. Suppose that RP ∋ θ 7→ ψθ = ηφθ (with η as in Theorem 3.2) is a parametrized family of compactly supported C 2 functions and that, for µ ∈ RP , 0 is a regular value of φµ . Then, one can establish via first order Taylor expansion that E[|Zµ |p ] < ∞ for every p < 32 , which is consistent with the lower bound from Theorem 2.4. Since by Thom’s transversality theorem [21, Chapter 3, Theorem 2.1] a generic C 2 function has 0 as a regular value, it follows – roughly speaking – that generically the variable Zθ has p moments for every p < 23 but as soon as ψν has higher order zeroes, the moment threshold drops to p < 45 . In other words, the generic critical threshold is 23 but the worst-case critical threshold is 54 . Remark 1.7. Using Malgrange’s preparation theorem [28], one can extend Theorem 3.2 to the case where the family φθ (x) is merely C ∞ and not real-analytic, as long as the functions φθ do not posses ”flat zeros” where all derivatives vanish. Remark 1.8. The assumption of Θ being compact in Theorem 3.2 is needed to achieve a uniform bound of E[|Zθ |p ] over all parameters θ ∈ Θ. However, the optimization algorithms presented and analyzed in Section 4 assume unconstrained optimization over θ ∈ RP and therefore – strictly speaking – the families of Theorem 3.2 do not fall into the unconstrained optimization framework studied in Section 4. This is not really an issue, since one can easily extend the families (ψθ )θ∈Θ to ψ̃θ := ψg(θ) , θ ∈ RP and g : RP → Θ smooth with uniformly bounded derivatives. It is easy to see that this new family still satisfies the same moment assumptions and, in particular, the Assumptions 4.2 needed in the convergence analysis of Section 4. Remark 1.9. The result of Theorem 3.2 requires the functions ψθ to be compactly supported, which is a restriction. However, the restriction is quite mild for applications related to Variational Monte Carlo, since it is known that bound states (e.g., eigenfunctions of H corresponding to eigenvalues below the ionization threshold) are exponentially decaying, see [17, Theorem 8.6] or [2]. Extending Theorem 3.2 to globally supported functions seems challenging; a potential avenue could be techniques similar to those used in [16, Section 3]. We leave this to future work. Remark 1.10. Theorem 3.2 does not yet apply fully to families ψθ with cusps, which by Kato’s cusp conditions [22] always appear in ground states of the multi-electron Hamiltonian. We do however expect that with similar arguments one can extend our results and leave the incorporation of cusps to future work. 1.2.3

Robust Optimization

Having clarified the precise regime of moments for the gradient random variable, we continue with our third main contribution: a new robust variant of VMC that provably converges to a stationary point of L, both in expectation (see Theorem 4.5) and with high probability (see Theorem 4.11). This algorithm – coined PSClip-VMC – arises as a generalization of our recent work [29] where robust gradient estimators are constructed using a per sample clipping procedure. Roughly speaking, the resulting gradient estimators assumes n i.i.d. samples X1 , . . . , Xn ∼ Xθ and corresponding i.i.d. samples, W1 , . . . , Wn ∼ Wθ and E1 , . . . , En ∼ Eθ and approximates ∇L(θ) = 2E[Zθ ] by a clipped empirical mean ( ) 1 n 2X 2 X αk p γα,p,k (Ek Wk ) − 2 γβ,2,k (Ek ) γβ,2,k (Wj ) , where γα,p,k (v) := v · min 1, . n n −n |v| k=1

k̸=j

In Theorem 4.5 we show that under weak moment assumptions, consistent with the one proved in Section 3, using this robust estimator guarantees convergence in expectation to a stationary point. Specifically, we show p that for any small enough, constant step-size and any batch-size n ≥ M 2(p−1) the iterates of PS-Clip-VMC satisfy M 1 X E(|∇θ L(θm )|2 ) = O M −1 . M m=1 In Theorem 4.11 we then show that under the same assumptions, and for any δ ∈ (0, 1) the iterates satisfy 2(p−1) ! M 1 X (log(4M/δ) + 1/4) p 2 |∇θ L(θm )| = O . M m=1 M with probability 1 − δ. To the best of our knowledge, these results establish, for the first time, convergence of VMC to a stationary point under realistic moment assumptions. Remark 1.11. Clipping only the local energy observations is a common practice for stabilizing training in variational Monte Carlo [20, 19, 8, 10, 14, 12]. In Appendix B, we show that this is sufficient to guarantee convergence in expectation, but not convergence with high probability in the heavy-tailed regime. In practice, local energy clipping is usually combined with a global clipping, or norm constraint, of the stochastic gradient. It is possible that this combination is also sufficient to ensure convergence with high probability. However, the preliminary experiments in Section 5 indicate that incorporating per-sample clipping of the gradient random variables can improve the stability of practical VMC training. 1.2.4

Experiments

We complement our theoretical results with preliminary numerical experiments using FermiNet [31, 40] applied to atoms with up to 18 electrons. Our results in Section 4 suggest that, in the heavy-tailed regime of VMC, one should clip both the local energies and the individual log-gradient contributions before averaging. The experiments in Section 5 are intended to test this principle in a realistic setting. To this end, we train FermiNet on Sulfur and Argon, using the standard components of modern neural-network-based VMC such as MCMC sampling, KFAC updates, and centered local-energy clipping. In order to stay close to competitive practical VMC training, the algorithm used in the experiments is not a literal implementation of the estimator analyzed in Section 4. In particular, instead of clipping thresholds which increase with the sample index, we use constant centered clipping thresholds, in analogy with the clipping rule commonly used for the local energy. The principle of per-sample clipping however remains. Our experiments indicate that per-sample gradient clipping can improve the stability and performance of VMC. In particular, using only half the number of samples per batch typically employed to train FermiNet, the model trained with PS-Clip-VMC matches the best published energy for Argon and improves upon the best published energy for Sulfur; see Table 1. The robustness of PS-Clip-VMC is particularly evident in the training trajectory for Argon shown in Figure 1. The standard method exhibits a sharp increase in energy around step 55,000, and is not able to recover and converge to a lower energy after that. Atom

Standard

PS-Clip-VMC

Best Published

S Ar

-398.1089(3) -527.5304(4)

-398.1110(3) -527.5420(3)

-398.110b -527.5419a

Table 1: Final energies for Sulfur and Argon using a batch-size of 2048. The cited values are from a: DeepErwin [12], b: CAS + experimental corrections [7].

1.3

Related Work

The literature on VMC is quite vast, see [5] and references therein. Especially for the Born-Oppenheimer multi-electron Hamiltonian the method has found renewed interest due to the empirical success of FermiNet 5

Figure 1: Training trajectories of Sulfur and Argon using a batchsize of 2048. The values are smoothed using a rolling mean of the last 2000 iterations. and related methods, see [31, 19, 34, 37, 35, 36, 38, 12, 13, 20, 8, 9, 10] and the survey paper [11]. The mathematical study of these methods is still at its infancy. Important results concern the expressivity of Deep Learning VMC ansätze, see for example [30, 27]. The optimization problem for VMC has been studied mathematically in [1, 26], albeit under moment conditions which – as we show – are not even satisfied for generic standard Slater-Jastrow ansätze. Regarding the existence of heavy tails, the most closely related work we found is the paper [44], where asymptotic expansions of the wavefunction and local energy are provided. Similarly to our first main result, Theorem 2.4, it is concluded that the energy random variable Eθ is heavy-tailed. Therefore, to our knowledge, [44] should be credited with first discovering the heavy-tailed nature of VMC algorithms. We note, however, that the main focus of [44] is to establish CLT-like results in order to construct asymptotic confidence intervals. Moreover, we emphasize that the results in [44] are non-rigorous and do not prove the validity of its asymptotic expansions (which would indeed require complicated arguments along the lines of our Lemma 2.11). Finally, unlike our work, [44] does not establish the existence of a certain number of moments or any rigorous convergence theorems.

1.4

Notation

For K ∈ N we write [K] := {1, . . . , K}. Moreover for p ∈ [1, ∞], k, d ∈ N and Ω ⊂ Rd with smooth boundary we write Lp (Ω), H k (Ω), C k (Ω), and Lploc (Ω) for the Lebesque space with integrability p, the k-th order Sobolev space, the space of k-times continuously differentiable functions, and the space of locally p-integrable functions. For N, d ∈ N and x = (x1 , . . . , xN ) ∈ Rd×N we occasionally write x̂i := (x1 , . . . , xi−1 , xi+1 , . . . , xN ) ∈ Rd×(N −1) .

1.5

Outline

The outline is as follows. In Section 2 we prove the non-existence of moments for general Slater-Jastrow ansätzte. In Section 3 we prove the existence of a minimal number of moments for real-analytic ansätze. In Section 4 we introduce and analyze novel optimization algorithms and prove their convergence. Finally, in Section 5 we present computational experiments.

Non-existence of Moments for a Slater-Jastrow Ansatz

The present section is devoted to proving that the random variables Eθ for the energy and Zθ for the energy gradient are generically heavy tailed.

2.1

Formulation of the Main Result

We start with describing the main settings of this section, starting with our main notation and assumptions for hamiltonians and associated random quantities. 6

Setting 2.1 (Hamiltonian). For d ∈ N, Ω ⊂ Rd open and ψ ∈ H 2 (Ω) let V ∈ L1loc (Ω) be such that the Hamiltonian 2 H (Ω) → L2 (Ω) H: ψ 7→ −∆ψ + Vψ is bounded (see for example [23, Theorem 4.3] for suitable conditions for this to hold). |ψ(x)|2 For x ∈ Rd and ψ ∈ H 2 (Ω) we define the local energy EL,ψ (x) := Hψ(x) ψ(x) and the density pψ (x) := ∥ψ∥2 . L2

For Xψ ∼ pψ we call the random variable Eψ := EL,ψ (Xψ ) the energy random variable. It satisfies that E[Eψ ] = ∥ψ∥12 ⟨Hψ, ψ⟩L2 . L2

If RP ∋ θ 7→ ψθ ∈ H 2 (Ω) is a parametrized model such that the map θ 7→ ψθ ∈ L2 (Ω) is differentiable, θ (Xθ ) we denote pθ := pψθ , Xθ := Xψθ , EL,θ := EL,ψθ , Eθ := Eψθ , and Wθ := ∇ψθ θψ(X . We call the random θ) variable Zθ := Eθ · Wθ − E[Eθ ] · E [Wθ ] the energy gradient random variable. It satisfies that 2E[Zθ ] = ∇θ E[Eθ ]

and

E [Wθ ] =

1 ∇θ ∥ψθ ∥2L2 , 2 ∥ψθ ∥2L2

see [4, Section 2.2]. The present section focuses on fermionic many-particle systems which are formally introduced in the following setting. Setting 2.2 (Fermionic Many-Particle Systems and Slater Determinants). Assume Setting 2.1. For N ∈ N, A ⊂ R3 and F (AN ) a function space on AN (such as L2 (AN ) or C(AN )) we denote n o Fa AN := ψ ∈ F AN : ∀σ ∈ SN , x = (x1 , . . . , xN ) ∈ AN : ψ(xσ(1) , . . . , xσ(N ) ) = (−1)sgn(σ) ψ(x1 , . . . , xN ) , where SN denotes the set or permutations of [N ] and for a permutation σ ∈ SN we denote with sgn(σ) its signature. If in Setting 2.1 it holds that Ω = AN and H is restricted to Ha2 (Ω), we speak of a fermionic many-particle system. For h1 , . . . , hN ∈ C(R3 ) we define the Slater determinant |h1 , . . . , hN ⟩ ∈ Ca (R3×N ) via   h1 (x1 ) . . . hN (x1 )   .. 3×N .. |h1 , . . . , hN ⟩(x1 , . . . , xN ) := det  ... . (5)  , (x1 , . . . , xN ) ∈ R . . h1 (xN )

...

hN (xN )

For M ∈ N and R = (R1 , . . . , RM ) ∈ R3×M define the collision set Σ(R) := {(x1 , . . . , xN ) ∈ R3×N : ∃i ∈ [N ], j ∈ [N ] \ {i} : xi = xj } ∪

[

(R3 )i−1 ×

i∈[N ]

[

{Rj } × (R3 )N −i

j∈[M ]

and assume that V ∈ C(R3×N \ Σ(R)). Finally, we introduce the Slater-Jastrow ansatz for the numerical solution of (2). Setting 2.3 (Slater-Jastrow Ansatz with variable exponent Slater-Type Orbitals (STOs)). Assume Settings 2.1 and 2.2. For l ∈ N0 let Yl,−l , . . . , Yl,l be spherical harmonics, e.g., an ONB of Hl := {u ∈ L2 (S2 ) : u harmonic homogeneous polynomial of degree l } ⊂ L2 (S2 ). Given L, I, K ∈ N, α1 , . . . , αI ∈ [0, ∞) the corresponding Slater-type orbitals (STOs) with variable exponents (see for example [18, Section 6.5.6]) are defined as x R3 ∋ x 7→ hαi ,l,m,k (x) := Yl,m |x|l+k e−αi |x| , l = 0, . . . , L, m = −l, . . . , l, k = 0, . . . , K−1, i = 1, . . . , I. |x| (6) 7

For α ∈ [0, ∞)I , c = (cl,m,k,i )l=0,...,L, m=−l,...,l, k=0,...,K−1, i=1,...,I ∈ R(L+1) KI and x ∈ R3 we denote h(α,c) (x) :=

cl,m,k,i hαi ,l,m,k (x).

l,m,k,i

Let Σ := Σ(0) and B, J, N ∈ N and RJ ∋ β 7→ Jβ ∈ C 2 (R3×N \ Σ) smooth with Jβ symmetric and N B 2 2 ∼ positive. Then, for ((αba , cba ))a=1,...,N, b=1,...,B ∈ [0, ∞)I × R(L+1) KL = [0, ∞)N BI × RN B(L+1) KL 2 and θ := β, ((αba , cba ))a=1,...,N, b=1,...,B ∈ RJ × [0, ∞)N BI × RN B(L+1) KL we denote ! B X ψθ (x) := Jβ (x) · |h(αb1 ,cb1 ) , . . . , h(αbN ,cbN ) ⟩(x) b=1

as the Slater-Jastrow ansatz with variable exponent Slater-type orbitals. We can now rigorously formulate our first main theorem. Theorem 2.4. Assume Settings 2.1 and 2.3. Assume further that L ≥ 1 and KI ≥ 3N − 3. Then, the set n h i o 2 θ ∈ RJ × [0, ∞)N BI × RN B(L+1) KL : E |Eθ |3 = ∞ ∧ E |Zθ |3/2 = ∞ (7) 2

contains a dense open subset in RJ × [0, ∞)N BI × RN B(L+1) KL . In other words, nonexistence of a 3rd moment of Eθ and of a 3/2-th moment of Zθ constitutes a generic property.

2.2

Proof of the Main Result

We now turn to the proof of Theorem 2.4. At a high level we will show that every ψθ has an arbitrarily small perturbation with certain regular zeros. Then, the existence of such zeros implies a lack of moments, as the following elementary results show. Lemma 2.5. Assume Setting 2.1, let A ⊂ Ω open with V ∈ C(A) and ψ ∈ C 2 (A). Assume there exists x ∈ A so that ψ(x) = 0, ∇ψ(x) ̸= 0, and ∆ψ(x) ̸= 0. (8) Then, it holds that E |Eψ |3 = ∞. If ψ = ψθ and additionally ∇θ ψθ (x) ̸= 0 with ∇θ ψθ continuous at x, then it also holds that h i E |Zθ |3/2 = ∞ .

(9)

(10)

Proof. Let d = 3n and x ∈ A with ψ(x) = 0 and ∇ψ(x) ̸= 0 and ∆ψ(x) ̸= 0. For y ∈ Rd we denote yd−1 := (y1 , . . . , yd−1 ). Without loss of generality, suppose that ∂x∂ d ψ(x) ̸= 0 and ∥ψ∥L2 (Ω) = 1. Then, by the implicit function theorem there exists an open set U ⊂ Rd−1 with xd−1 ∈ U and g ∈ C 1 (U ) with g(xd−1 ) = xd

∧

∀ yd−1 ∈ U : ψ(yd−1 , g(yd−1 )) = 0.

(11)

For ε > 0 denote Bε := {(yd−1 , g(yd−1 ) + τ ) : yd−1 ∈ U, τ ∈ (−ε, ε)} ⊂ Rd . Then, Bε is open and for every integrable φ : Bε → R it holds that Z Z εZ φ(y)dy = φ(yd−1 , g(yd−1 ) + τ )dyd−1 dτ. Bε

−ε

By noting that V(x)ψ(x) = 0 and shrinking U and ε if necessary, we can ensure that Bε ⊂ A and inf |Hψ(y)| =: δ > 0.

y∈Bε

(12)

Now, note that 3

Hψ(y) |ψ(y)|2 dy ψ(y) Ω Z 3 ≥δ |ψ(y)|−1 dy Bε Z εZ |ψ(yd−1 , g(yd−1 ) + τ )|−1 dyd−1 dτ. = δ3

E |Eψ |3 =

−ε

(13) (14) (15)

Furthermore, using (11), it holds for yd−1 ∈ U and τ ∈ (−ε, ε) that ψ(yd−1 , g(yd−1 ) + τ ) = ψ(yd−1 , g(yd−1 )) + τ · =τ·

∂ ψ(yd−1 , g(yd−1 )) + η(yd−1 , τ ) ∂xd

with |η(yd−1 , τ )| ≤ τ 2 ∥ψ∥C 2 (Bε ) . Therefore, for all yd−1 ∈ U and τ ∈ (−ε, ε) we have that |ψ(yd−1 , g(yd−1 ) + τ )| ≤ |τ |(1 + ε)∥ψ∥C 2 (Bε ) . Inserting this inequality into (15) yields that E |Eψ |3 ≥

|U |δ 3 (1 + ε)∥ψ∥C 2 (Bε )

Z ε

1 dτ = ∞, −ε |τ |

which proves (9). To prove (10), we can apply the same reasoning to the integral Z Bε

Hψ(y) |∇θ ψ(y)| ψ(y) ψ(y)

3/2

|ψ(y)|2 dy ,

(16)

and use the fact that (by the assumption that ∇θ ψθ (x) ̸= 0 and shrinking U, ε if needed) it holds that |∇θ ψ(y)| > 0 for all y ∈ Bε . This implies that Z Bε

Hψ(y) |∇θ ψ(y)| ψ(y) ψ(y)

3/2

|ψ(y)|2 dy = ∞.

(17)

Since (17) is the dominant term in Zθ , it follows that E[|Zθ |3/2 ] = ∞. Lemma 2.6. Let A ⊂ R3×N open. Then, O := {ψ ∈ C 2 (A) : ∃x ∈ A : ψ(x) = 0, ∇ψ(x) ̸= 0, and ∆ψ(x) ̸= 0} is open in the local (compact-open) topology of C 2 (A). If ψ = ψθ such that the mapping RP ∋ θ 7→ (ψθ , ∇θ ψθ ) ∈ C 2 (A) × C(A, RP ) is continuous in the local (compact-open) topology, then the set O′ := {θ ∈ RP : ∃x ∈ A : ψθ (x) = 0, ∇ψθ (x) ̸= 0, ∆ψθ (x) ̸= 0, and ∇θ ψθ (x) ̸= 0} is open. Proof. Recall that in the local C 2 topology, for every ψ ∈ C 2 (A), every compact B ⊂ A and for every ε > 0, the set V (B, ε, ψ) := {φ ∈ C 2 (A) : ∥φ − ψ∥C 2 (B) < ε} is open. We will show that for every ψ ∈ O there is B, ε such that V (B, ε, ψ) ⊆ O. 9

Let ψ ∈ O with x as in the definition of O. Let B ⊂ A be a τ -neighborhood of x and let y ∈ R3×n be a unit vector with s := y · ∇ψ(x) > 0. Then, it holds for all t ∈ [−τ, τ ] that φ(x + ty) = φ(x) + ts + ηt , where |ηt | ≤ ∥φ∥C 2 (B) · t2 . By choosing τ<

s ∥ψ∥C 2 (B)

and

ε<

τ s − ∥ψ∥C 2 (B) τ 2 1 + τ + τ2

it is easy to see that, whenever φ ∈ V (B, ε, ψ), we have that φ(x − τ y) < 0

and

φ(x + τ y) > 0

and therefore, there is τ∗ ∈ (−τ, τ ) such that with x∗ := x + τ∗ y ∈ B it holds that φ(x∗ ) = 0. By reducing τ, ε further if needed and using that ∇ψ(x) ̸= 0 and ∆ψ(x) ̸= 0 we can ensure that ∇φ(x∗ ) ̸= 0

and

∆φ(x∗ ) ̸= 0.

This proves that V (B, ε, ψ) ⊂ O as needed. The statment on openness of O′ follows now immediately from the continuity of the mapping θ 7→ (ψθ , ∇θ ψθ ). The key difficulty in the proof of Theorem 2.4 is now to show that any ψθ can be slightly perturbed so that a zero as in Lemma 2.5 exists. We will ultimately achieve this by slightly changing the underlying Slater-type orbital, but our results will require sufficient expressivity of the underlying Slater-type orbitals. A key condition in this context is the following. Definition 2.7. Let Vr ⊂ C 2 (0, ∞) be such that for m ∈ N, every 0 < τ1 < · · · < τm and every (i) (i) di (νj )j∈[m], i∈{0,1,2} there is v ∈ Vr with dt for all j ∈ [m], i ∈ {0, 1, 2}. Then, Vr is said i v(τj ) = νj to satisfy the m-Hermite interpolation property (m-HIP). We also mention that the m-HIP is closely related to Vr constituting an extended Chebychev system [32] – a property that is known to be satisfied by exponentially weighted polynomials, as the following lemma shows. We emphasize however, that the m-HIP only applies to the radial part of a 3-dimensional orbital. Lemma 2.8. If Vr = span{rk eαi r : i = 1, . . . I, k = 0, . . . , mi − 1} with αi ∈ R pairwise disjoint and PI i=1 mi ≥ 3m, then Vr satisfies the m-HIP. m 1 m I d d Proof. Observe that Vr = ker dt − α1 . . . dt − αI . By [3, Theorem 9], the space Vr is an extended Chebychev system in the sense of [32, Definition III.1.8], which implies the m-HIP. The following Lemma establishes that the (N − 1)-HIP implies, for specific nodes xi ∈ R3 , i ∈ [N ], that a multiplicative linear factor allows the construction of 3-dimensional orbitals with prescribed point values and derivatives. We mention that such a result is by no means obvious: Due to the well-known Mairhuber-Curtis Theorem [45], general interpolation properties of linear subspaces for d-dimensional functions are in general highly nontrivial for d ≥ 2. Lemma 2.9. Let N ∈ N, suppose that Vr ⊂ C 2 (0, ∞) satisfies the (N − 1)-HIP and let V := {R3 ∋ x 7→ (w · x)vr (|x|) : w ∈ R3 ∧ vr ∈ Vr }. Then, for every pairwise distinct x1 , . . . , xN ∈ R3 \ {0} with x1 , x2 non-collinear such that x1,1 ̸= 0

∧

|x1 | = |x2 |

∧

∀i ∈ {3, . . . , N }, j ∈ [N ] \ {i} : |xi | ̸= |xj | ,

and for every ν (1) , ν (2) ∈ R there exists v ∈ V with 10

(18)

(i) v(xi ) = 0 for all i ∈ [N ], (1)

(ii) ∇v(x1 ) = xν 1,1 x1 and ∇v(xi ) = 0 for i = {2, . . . , N }, and (iii) ∆v(x1 ) = ν (2) and ∆v(xi ) = 0 for i = {2, . . . , N }. Proof. Observe that for any v(x) = (w · x)vr (|x|) with vr ∈ Vr and k ∈ {1, 2, 3}, it holds that xk 4 ′ ′ ′′ ∂k v(x) = wk vr (|x|) + (w · x)vr (|x|) and ∆v(x) = (w · x) vr (|x|) + v (|x|) . |x| |x| r

(19)

Now, pick vr ∈ Vr and w ∈ R3 so that ∀i ∈ [N ] : vr (|xi |) = 0 ∧ ∀i ∈ {3, . . . , N } : vr′ (|xi |) = vr′′ (|xi |) = 0 ∧ w · x1 =: τ ̸= 0 ∧ w · x2 = 0. (20) Note that this is possible due to the (N − 1)-HIP, and we are still free to choose values for vr′ (|x1 |), vr′′ (|x1 |). Using (19) and (20), it is easy to check that v(x) := (w · x)vr (|x|) satisfies ∀i ∈ [N ] : v(xi ) = 0

∧

∀i ∈ {2, . . . , N } : ∆v(xi ) = 0 ∧ ∇v(xi ) = 0.

Furthermore, it holds that x1,k ∂k v(x1 ) = τ vr′ (|x1 |) |x1 |

and

∆v(x1 ) = τ

vr′′ (|x1 |) +

4 ′ v (|x1 |) . |x1 | r

By choosing the values vr′ (|x1 |), vr′′ (|x1 |) appropriately, we can ensure that (i), (ii) and (iii) are satisfied. The following key lemma constructs interpolation nodes which satisfy three properties simultaneously: the assumptions of Lemma 2.9, being a zero of a given antisymmetric function ψ, as well as a non-vanishing property of an associated Slater determinant. The proof uses a topological argument similar to the one in the proof of the Mairhuber-Curtis theorem. Lemma 2.10. Assume Setting 2.2 and let v2 , . . . , vN ∈ C(R3 ) linearly independent and ψ ∈ Ca (R3×N ). Then, for R = (R1 , . . . , Rm ) ∈ R3×M there exist x1 , . . . , xN ∈ R3×N \ Σ(R) with x1 , x2 non-collinear and x1,1 ̸= 0

∧

|x1 | = |x2 |

∧

∀i ∈ {3, . . . , N }, j ∈ [N ] \ {i} : |xi | ̸= |xj |

(21)

such that (i) ψ(x1 , . . . , xN ) = 0, and (ii) |v2 , . . . , vN ⟩(x2 , . . . , xN ) ̸= 0. Proof. Since v2 , . . . , vN are linearly independent, there exist z = (z∗ , z3 , . . . , zN ) ∈ R3×(N −1) so that |v2 , . . . , vN ⟩(z∗ , z3 , . . . , zN ) ̸= 0. Due to continuity, there exists an open neighborhood U of z in R3×(N −1) so that |v1 , . . . , vN ⟩ vanishes nowhere on U . In particular, we can perturb z so that |zi | ̸= |zj |

∀i ∈ {∗, 3, . . . , N }, j ∈ {∗, 3, . . . , N } \ {i} and |zi | ∈ /

M [

{|Rl |}

∀i ∈ {∗, 3, . . . , N },

l=1

as well as z∗,1 ̸= 0. Now, pick a circle C = {γ(t) = z∗ + ε sin(t)a + ε cos(t)b : t ∈ [0, 2π)} with unit vectors a, b ∈ {z∗ }⊥ and C × {z3 } × . . . {zN } ⊂ U . By making ε sufficiently small, and using that z∗,1 ̸= 0, we can ensure that ∀x ∈ C : x1 ̸= 0 and ∀x, y ∈ C : (x ̸= y) ⇒ x, y not collinear. Note that ∀t ∈ [0, 2π) : |γ(t)|2 = |z∗ |2 + ε2 . 11

(22)

Now, consider the function τ (t) := ψ(γ(t), γ(t + π), z3 , . . . , zN ⟩ Due to antisymmetry of ψ, it holds that τ (0) = −τ (π). Since τ is continuous, there must exist t0 ∈ [0, π) such that ψ(γ(t0 ), γ(t0 + π), z3 , . . . , zN ) = 0. In summary, (x1 , x2 , x3 , . . . , xN ) = (γ(t0 ), γ(t0 + π), z3 , . . . , zN ) satisfies the claimed properties. The following lemma constitutes the main technical result of this section and establishes the existence of perturbations with regular zeros. Lemma 2.11. Assume Settings 2.2 and 2.3. Suppose that Vr ⊂ C 2 (0, ∞) satisfies the (N − 1)-HIP and let V be a vector space with V ⊇ {R3 ∋ x 7→ (w · x)vr (|x|) : w ∈ R3 ∧ vr ∈ Vr }. Then, for any linearly independent functions v1 , . . . , vN ∈ C(R3 ) with v1 ∈ V , any F ∈ Ca2 (R3×N \ Σ), any symmetric J ∈ C 2 (R3×N \ Σ) with J > 0, there exist x ∈ R3×N \ Σ, vδ , u ∈ V \ {0} and λ0 > 0 such that for all λ ∈ (0, λ0 ) and ψδ,λ := F + J · |v1 + λvδ , v2 , . . . , vN ⟩ we have that ψδ,λ (x) = 0,

∇ψδ,λ (x) ̸= 0,

and

∆ψδ,λ (x) ̸= 0

(23)

and 1 (F (x) + J (x) · |v1 + λvδ + hu, v2 , . . . , vn ⟩(x) − ψδ,λ (x)) = J (x) · |u, v2 , . . . , vn ⟩(x) ̸= 0. h→0 h

(24)

lim

Proof. Write ψ := F + J · |v1 , . . . , vN ⟩ and ψδ := ψδ,1 . For i ∈ [N ] and (x1 , . . . , xN ) ∈ R3×N let x̂i := (x1 , . . . , xi−1 , xi+1 , . . . , xN ) ∈ R3×(N −1) . Then, for all x = (x1 , . . . , xN ) ∈ R3×N it holds that ψδ (x) = ψ(x) + J (x)

N X

(−1)i+1 vδ (xi )|v2 , . . . , vN ⟩(x̂i )

(25)

i=1 N X ∂ ∂ ∂ (−1)i+1 vδ (xi ) ψδ (x) = ψ(x) + J (x) ∂1 vδ (x1 )|v2 , . . . , vN ⟩(x̂1 ) + |v2 , . . . , vN ⟩(x̂i ) ∂x1,1 ∂x1,1 ∂x1,1 i=2

N X ∂ (−1)i+1 vδ (xi )|v2 , . . . , vN ⟩(x̂i ) (26) + J (x) ∂x1,1 i=1   X ∆i δψ(x) = ∆i ψ(x) + J (x)  (−1)j+1 vδ (xj )∆i |v2 , . . . , vN ⟩(x̂j ) + (−1)i+1 ∆vδ (xi )|v2 , . . . , vN ⟩(x̂i ) j̸=i

+ ∆i J (x)

N X

! i+1

(−1)

vδ (xi )|v2 , . . . , vN ⟩(x̂i )

i=1





+ 2∇i J (x) · (−1)i+1 ∇vδ (xi )|v2 , . . . , vN ⟩(x̂i ) +

(−1)j+1 vδ (xj )∇i |v2 , . . . , vN ⟩(x̂j ) . (27)

j̸=i

Now, pick x = (x1 , . . . , xN ) according to Lemma 2.10 and, using these points, vδ according to Lemma 2.9 with ν (1) , ν (2) to be determined. Then, with τ := |v2 , . . . , vN ⟩(x̂1 ) ̸= 0, J (x) ̸= 0 and noting that (1) ∇vδ (x1 ) = νx1,1x1 the formulas (25), (26), (27) simplify to ψδ (x) = 0

(28)

∂ ∂ ψδ (x) = ψ(x) + ν (1) J (x)τ ∂x1,1 ∂x1,1 ∆ψδ (x) = ∆ψ(x) + ν (2) J (x)τ +

2ν (1) (∇1 J (x) · x1 ) . x1,1

(29) (30)

Since we can choose ν (1) , ν (2) freely, it is easy to ensure that ∂x∂1,1 ψδ (x) ̸= 0 and ∆ψδ (x) ̸= 0. By scaling vδ (and hence (ν (1) , ν (2) )), we can ensure via (29), (30) that ∂x∂1,1 ψδ,λ (x) ̸= 0 and ∆δψλ (x) ̸= 0, whenever λ > 0 is sufficiently small. This proves (23). It remains to prove (24). To this end, let u ∈ V with u(x1 ) ̸= 0 and u(xi ) = 0 for all i = 2, . . . , N (such u can be constructed by picking vr to vanish at |x3 |, . . . , |xN | but not on |x1 | = |x2 |, and w ∈ R3 with w · x2 = 0 and w · x1 ̸= 0. Then, u(x) := (w · x)vr (|x|) satisfies the desired properties). Hence, it holds that |u, v2 , . . . , vN ⟩(x) =

N X

(−1)1+i u(xi )|v2 , . . . , vN ⟩(x̂i ) = u(x1 )τ ̸= 0.

i=1

Since J (x) ̸= 0 and τ ̸= 0, this proves the inequality in (24). The equality in (24) is immediate from the multilinearity of the determinant. Using Lemma 2.11, we are finally in a position to prove Theorem 2.4. Proof of Theorem 2.4. Since the map θ 7→ (ψθ , ∇θ ψθ ) is continuous with respect to the local topology, by Lemma 2.6 the set 2

O′ := {θ ∈ RJ ×[0, ∞)N BI ×RN B(L+1) KL : ∃x ∈ R3×N : ψθ (x) = 0, ∇ψθ (x) ̸= 0, ∆ψθ (x) ̸= 0, and ∇θ ψθ (x) ̸= 0} 2

is open in RJ × [0, ∞)N BI × RN B(L+1) KL . We will show that O′ is also dense in RJ × [0, ∞)N BI × 2 RN B(L+1) KL . To this end, let ! B X ψθ = Jβ · |h(αb1 ,cb1 ) , . . . , h(αbN ,cbN ) ⟩ . b=1 B are linearly By removing a nowhere dense set of parameters we may assume that the functions vi := hαB i ,ci 3 3×N B N independent in C(R ) (to see this we note that for (x1 , . . . , xN ) ∈ R the mapping Φ : (αB i , ci )i=1 7→ B B N det (vi (xj )) is real analytic and not identically zero, therefore the set of parameters (αi , ci )i=1 with B N B I Φ (αB i , ci )i=1 = 0 is nowhere dense) and that α1 = (α1 , . . . , αI ) ∈ [0, ∞) has pairwise distinct entries. Consider the finite dimensional vector space

V := span{hαB : l = 0, . . . , L, m = −l, . . . , l, k = 0, . . . , K − 1, i = 1, . . . , I} 1,i ,l,m,j and observe that, since h1 = {R3 ∋ x 7→ w · x : w ∈ R3 } and L ≥ 1 it holds that V ⊇ {R3 ∋ x 7→ (w · x)vr (|x|) : w ∈ R3 ∧ vr ∈ Vr }, where

Vr = span{rk e−α1,i r : i = 1, . . . I, k = 0, . . . , K − 1}. Since by assumption KI ≥ 3(N − 1), it follows from Lemma 2.8 that Vr satisfies the (N − 1)-HIP. Therefore, B B in the parameter c , we can apply Lemma 2.11 with J = Jβ and F = using the linearity of hαB 1 1 ,c1 PB−1 Gβ · b=1 |hαb1 ,cb1 , . . . , hαbN ,cbN ⟩, to conclude that an arbitrarily small perturbation θ′ of θ produces x ∈ R3×N with ψθ′ (x) = 0, ∇ψθ′ (x) ̸= 0 and ∆ψθ′ (x) ̸= 0. (31) Moreover, by (24) there exists u ∈ V with Jβ (x) · |u, v2 , . . . , vN ⟩(x) ̸= 0. By noting that ∇cB ψθ′ (x) ̸= 0 ⇔ ∃u ∈ V : |u, v2 , . . . , vN ⟩(x) ̸= 0 1 we conclude that also ∇θ ψθ′ (x) ̸= 0 ′

′

(32) ′

holds true. Together, (31) and (32) establish that θ ∈ O and therefore O is dense. Finally, by Lemma 2.5 it holds that n h i o 2 O′ ⊂ θ ∈ RJ × [0, ∞)N BI × RN B(L+1) KL : E |Eθ |3 = ∞ ∧ E |Zθ |3/2 = ∞ . The theorem is proven. 13

Existence of Moments

In this section, we establish the key fact that, whenever the functions ψθ are real-analytic, a guaranteed minimal number of moments exists for the random variables Eθ and Zθ .

3.1

Formulation of the Main Results

Before stating and proving our main result, we begin by noting the following trivial moment properties. They will be needed in the analysis of the robust optimization methods introduced in Section 4. Theorem 3.1. Assume Setting 2.1. Then, for every ψ ∈ H 2 (Ω) \ {0} it holds that ∥Hψ∥2L2 < ∞. E |Eψ |2 = ∥ψ∥2L2

(33)

For θ 7→ ψθ ∈ H 2 (Ω) \ {0} with θ 7→ ψθ differentiable in L2 (Ω) it holds that ∥∇θ ψθ ∥2L2 E |Wθ |2 = < ∞. ∥ψ∥2L2

(34)

Proof. This follows directly from the definition. We now state the main result of this section. Theorem 3.2. Assume Setting 2.1. Let K ⊂ Ω compact with smooth boundary, ε > 0, P ∈ N and P ⊂ RP open so that the mapping P × Ω ∋ (θ, x) 7→ φθ (x) is in C ω (P × Ω) and Q ⊂ P compact. Assume further that 5/4 V ∈ Lloc . Then, for any η ∈ C ∞ (Ω) with supp(η) = K it holds for the parametrized family θ 7→ ψθ := ηφθ that 5 (35) sup{p ∈ (0, ∞) : sup E[|Zθ |p ] < ∞} ≥ . 4 θ∈Q Remark 3.3. A function η ∈ C ∞ (Ω) satisfying the properties in the statement of Theorem 3.2 can for example be constructed as follows. First, note that for ε > 0 sufficiently small and for τ > 0 defining Aτ := {x ∈ Rd : dist(x, ∂K) ≤ τ } ⊂ Rd the distance function dist(x, ∂K) is in C ∞(Aε/2 ∩ K). Now, let 1 s ∈ C ∞ (K) with s|K\Aε ≡ 1 and s|K∩Aε/2 = dist(·, ∂K) and let η(x) := exp(1) · exp − s(x) · 1K (x). It is

a standard fact that η ∈ C ∞ (Ω) (this can be seen from the fact that the function t 7→ e−1/t · 1(0,∞) is in C ∞ (R)). Furthermore, it holds that supp(η) = K and ∀x ∈ K : dist(x, ∂K) > ε ⇒ η(x) = 1. We finally note that the energy gradient of real-analytic parametrizations are Lipschitz continuous. This property will be needed to prove convergence of the robust optimization methods introduced in Section 4. Theorem 3.4. Assume Setting 2.1 and consider the mapping θ 7→ ψθ from Theorem 3.2. Assume additionally that ψθ ̸= 0 for every θ ∈ Q. Then, the mapping θ 7→ ψθ is globally Lipschitz continuous as a map from Q to H 2 (Ω) and continuously differentiable as a map from Q to L2 (Ω), with globally Lipschitz continuous derivative. Moreover, there exists a constant C ∈ (0, ∞) such that ∀θ1 , θ2 ∈ Q :

|E[Zθ1 ] − E[Zθ2 ]| ≤ C|θ1 − θ2 |.

(36)

Proof. We first show (using the notation of Theorem 3.2) that θ 7→ ψθ is globally Lipschitz continuous as a map from Q to H 2 (K) ⊂ H 2 (Ω). Let f (θ, x) := φθ (x). By assumption, this is a real-analytic function and consequently all its derivative, including ∇f, ∇2 f, ∇3 f , are also real-analytic on P × Ω (and therefore bounded on compact sets). Since Q × K is compact this directly implies that there exists a constant C ∈ (0, ∞) with ∥f (θ1 , ·)−f (θ2 , ·)∥L∞ (K) +∥∇f (θ1 , ·)−∇f (θ2 , ·)∥L∞ (K) +∥∇2 f (θ1 , ·)−∇2 f (θ2 , ·)∥L∞ (K) ≤ C|θ1 −θ2 |

∀θ1 , θ2 ∈ Q.

Since η ∈ C ∞ (K) this directly implies the existence of a (possibly different) constant C ∈ (0, ∞) with ∥ψθ1 − ψθ2 ∥H 2 (K) ≤ C|θ1 − θ2 | 14

∀θ1 , θ2 ∈ Q.

This shows that the map θ 7→ ψθ is globally Lipschitz as a map from Q to H 2 (K). Using similar arguments one can show that the mapping θ 7→ ψθ is infinitely often continuously differentiable as a map from Q to L2 (K), which in particular implies that there exists C ∈ (0, ∞) with ∥ψθ1 − ψθ2 ∥L2 (K) + ∥∇θ ψθ1 − ∇θ ψθ2 ∥L2 (K) ≤ C|θ1 − θ2 |

∀θ1 , θ2 ∈ Q.

It remains to show (36). To this end note that Z 1 1 ∆ψθ (x) · ∇θ ψθ (x) − ⟨Hψθ , ψθ ⟩ · ⟨∇θ ψθ , ψθ ⟩. E[Zθ ] = ∥ψθ ∥2L2 K ∥ψθ ∥4L2 By the previous considerations (and using that H : L2 → H 2 is bounded), this is Lipschitz continuous (in fact, it is easy to see that it is C ∞ ) in θ. Equation (36) follows.

3.2

Proof of the Main Result

We now prove Theorem 3.2. We will make extensive use of the well-known structure on the local geometry of nodal sets of analytic functions, as exhibited by the Weierstrass preparation theorem. Theorem 3.5 (Weierstrass Preparation Theorem). Let Ω ⊂ Rd open. Then, for each f ∈ C ω (Ω), and every x = (x1 , . . . , xd ) ∈ Ω such that ∀l = 0, . . . , k − 1 :

∂l f (x) = 0 ∂xl1

∧

∂k f (x) ̸= 0 ∂xk1

there is an interval I with x1 ∈ I and an open set V ⊂ Rd−1 with x̂1 ∈ V and I × V ⊂ Ω and unique functions ai ∈ C ω (V ), i = 0, . . . , k − 1, b ∈ C ω (I × V ) with ∀l = 0, . . . , k − 1 : al (x̂1 ) = 0 ∧ ∀y = (y1 , ŷ1 ) ∈ I × V : b(y) ̸= 0 ∧ f (y) =

(y1 − x1 ) +

k−1 X

! l

al (ŷ1 )(y1 − x1 )

· b(y).

(37)

l=0

Proof. A proof for holomorphic functions can be found in [25, Theorem 4.41], the real-analytic statement follows by noting that every real-analytic function can be locally extended to a holomorphic function. With the Weierstrass preparation theorem at hand, we can now reduce the local integrability Pk−1 of the random variables Eθ , Zθ to studying integrability properties of the Weierstrass polynomials y1 + l=0 al (yd−1 )y1l . This will involve quotients between gradients and Laplacians of the Weierstrass polynomial and the Weierstrass polynomial itself. Since applying the Laplacian requires us to differentiate the Weierstrass polynomial in every coordinate yi and since the root functions al are in general not analytic, we need to make sure that a Weierstrass preparation is possible in every coordinate. This is ensured by the following Lemma. Lemma 3.6. Let Ω ⊂ Rd open. Let f ∈ C ω (Ω), x = (x1 , . . . , xd ) ∈ Ω and ordf (x) := min{m ∈ N0 : ∃α := (α1 , . . . , αd ) ∈ Nd0 :

d X

αi = m ∧

i=1

∂ α1 ∂ αd f (x) ̸= 0}. α1 · · · d ∂xd ∂xα d

Then, there exist orthonormal vectors r1 , . . . , rd ∈ Rd such that ∀i ∈ [d] : ∀l = 0, . . . , ordf (x) − 1 :

dl f (x + tri )|t=0 = 0 dtl

∧

dordf (x) f (x + tri )|t=0 ̸= 0 . dtordf (x)

Proof. Since for l < ordf (x) all derivatives of f at x of order l vanish, it follows that for any r ∈ Rd and l < ordf (x) it holds that dl f (x + tr)|t=0 = 0. dtl 15

P ordf (x)! ∂ α1 ∂ αd α Furthermore, we note that with H(r) := α=(α1 ,...,αd )∈Nd : Pd αi =ordf (x) α! α ··· α f (x) · r , we ∂xd 1 ∂xd d 0 i=1 have that dordf (x) f (x + tr)|t=0 = H(r). (38) dtordf (x) Since this is a homogeneous polynomial, it follows that B := {r ∈ Sd−1 : H(r) = 0} O(d) → Sd−1 d−1 is closed and has empty interior in S . For i ∈ [d], let πi : and note that, R = (r1 , . . . , rd ) 7→ ri −1 since πi is continuous and open, the set πi (B) is closed and has empty interior in O(d). Therefore, there S ordf (x) exists R = (r1 , . . . , rd ) ∈ O(d) \ i∈[d] πi−1 (B) which by (38) satisfies that d ordf (x) f (x + tri )|t=0 ̸= 0 for all dt i ∈ [d]. Roughly speaking, the local integrability of Eθ , Zθ is determined by integrability properties of quantities related to Weierstrass polynomials. The following lemma establishes the necessary integrability properties for polynomials that will later be used in the proof of Theorem 3.2. Pk−1 Lemma 3.7. Let u(t) = tk + l=0 al tl , t ∈ I, and I ⊂ R a bounded interval. Then, for each p < 3/2 there is a constant C1 = C1 (I, maxk−1 l=0 |al |, p) with Z |u′ (t)|p |u(t)|2−2p dt ≤ C1 (39) I

and for each p < 5/4 there is a constant C2 = C2 (I, maxk−1 l=0 |al |, p) with Z |u′′ (t)|p |u(t)|2−2p dt ≤ C2 .

(40)

Proof. Consider the factorization u(t) =

k Y (t − αl ),

αl ∈ C.

l=1

Then, we have that ′

u (t) =

k X

u′′ (t) =

(t − αl ),

i=1 l∈[k]\{i}

(t − αl )

i∈[k], j∈[k]\{i} l∈[k]\{i,j}

and therefore k u′ (t) X = (t − αi )−1 , u(t) i=1

u′′ (t) = u(t)

(t − αi )−1 (t − αj )−1 .

i∈[k], j∈[k]\{i}

By the Cauchy root bound (together with the fact that the highest degree term in p has coefficient 1), it holds that k−1 |αl | ≤ 1 + max |al |. l=0

Hence, we obtain that u′ (t) u(t)

|u(t)|2−p ≤ k p−1

k X

|t−αi |2−2p

i=1

k−1

|t−αl |2−p ≤ k p−1 (|t|+1+ max |al |)(2−p)(k−1) l=0

l∈[k]\{i}

k X i=1

which is integrable on I if 2 − 2p > −1 or equivalently p < 23 . This proves (39). To establish (40) we note that with some constant Ck,p u′′ (t) u(t)

|u(t)|2−p ≤ Ck,p

|t − αi |2−2p |t − αj |2−2p

i∈[k] j∈[k]\{i}

Y l∈[k]\{i,j}

|t − αl |2−p .

|t−αi |2−2p ,

Hence, it is enough to bound terms of the form Z |t − αi |2−2p |t − αj |2−2p dt. I

By Cauchy-Schwarz, Z

2−2p

|t − αi |

2−2p

|t − αj |

Z dt ≤

4−4p

|t − αi |

1/2 Z 1/2 4−4p , dt |t − αj | dt

which is integrable on I if 4 − 4p > −1, or equivalently p < 45 . This, together with the estimate for |αl | above, proves (40). We are now ready to prove Theorem 3.2. Proof of Theorem 3.2. In view of the definition of Zθ we need to show that for p < 54 it holds that ∆ψµ (y) · ∇µ ψµ (y) ψµ (y)2 Ω

Z sup µ∈Q

|ψµ (y)|2 dy < ∞.

Denote f : (θ, x) 7→ φθ (x), and for (θ, x) ∈ Q × K fixed, let k = ordφθ (x). By Lemma 3.6, there are orthonormal vectors r1 , . . . , rd ∈ Rd with ∀i ∈ [d] : ∀l = 0, . . . , k − 1 :

dl f (θ, x + tri )|t=0 = 0 dtl

∧

dk f (θ, x + tri )|t=0 ̸= 0 . dtk

For ease of notation from now on we assume that r1 , . . . , rd is the standard Euclidean basis and consequently that ∂l ∂k ∀i ∈ [d] : ∀l = 0, . . . , k − 1 : f (θ, x) = 0 ∧ f (θ, x) ̸= 0. (41) ∂xli ∂xki Important Remark: Since ultimately we will develop estimates for ∆φθ and since ∆ is invariant under orthogonal coordinate transforms this reduction will be harmless. Let us introduce some convenient notation: For i ∈ [d] and y ∈ Rd we denote ŷi := (y1 , . . . , yi−1 , yi+1 , . . . , yd ) ∈ d−1 i−1 d−1 d−1 R and for I ⊂ R and a product set V = j=1 Vj ⊂ Rd−1 denote I ×i V := j=1 Vj × I × j=i+1 Vj . Due to (41), we can now apply Theorem 3.5 simultaneously in each coordinate of x (and, for ease of notation, assuming that x1 = 0), yielding for each i ∈ [d] an interval Ii with xi ∈ Ii , an open product set Vi ⊂ Rd−1 with x̂i ∈ Vi and Ii ×i Vi ⊂ Ω, an open set Wi ⊂ P with θ ∈ Wi , as well as unique functions ail ∈ C ω (Wi ×Vi ), l = 0, . . . , k − 1, bi ∈ C ω (Wi × Ii ×i Vi ) with

∀l = 0, . . . , k − 1 : ail (θ, x̂i ) = 0 i

∀(µ, y) ∈ Wi × Ii ×i Vi : b (µ, y) ̸= 0 ! k−1 X ∀(µ, y) ∈ Wi × Ii ×i Vi : f (µ, y) = yik + ail (µ, ŷi )yil · bi (µ, y).

(42) (43) (44)

l=0

By shrinking the sets Wi , Ii , Vi if necessary, we can in particular assume that |bi (µ, y)| > ϵ > 0 for all (µ, y) ∈ Wi × Ii ×i Vi and i ∈ [d]. We now aim to show that p

Z µ∈

sup T

d i=1 Wi

i=1 Ii ×i Vi

(45)

Since due to analyticity ∇µ f (µ, y) is locally bounded, (45) follows if we can establish that Z ∀i ∈ [d], l ∈ {1, 2} : sup µ∈Wi

Ii ×i Vi

∂l f (µ, y) ∂yil f (µ, y)2

|f (µ, y)|2 dy < ∞.

(46)

Pk−1 In order to prove (46), we use the representations (42), (43), and (44). Setting ui,µ,ŷi (yi ) := yik + l=0 ail (µ, ŷi )yil , it follows from (43) that (46) holds true if for all i ∈ [d]: Z Z p p u′i,µ,ŷi (yi ) u′′i,µ,ŷi (yi ) 2 sup |u (y )| |ui,µ,ŷi (yi )|2 dy < ∞ (47) dy < ∞ ∧ sup i i,µ,ŷ i 2 2 µ∈Wi Ii ×i Vi ui,µ,ŷi (yi ) µ∈Wi Ii ×i Vi ui,µ,ŷi (yi ) Finally, we note that, if p < 45 , (47) is now a direct consequence of Lemma 3.7 using the fact that the coefficient functions ail (µ, ŷi ) are bounded on Ii ×i Vi . We have thus established (45). Next, we note that Z p ∆ψµ (y) · ∇µ ψµ (y) |ψµ (y)|2 dy = sup Td 2 T ψ (y) µ I × V µ∈ d W i i i i i=1 i=1 Z p ∆(η(y)f (µ, y)) · η(y)∇µ f (µ, y) sup |η(y)f (µ, y)|2 dy < ∞ (48) Td 2 f (µ, y)2 T η(y) I × V µ∈ d Wi i=1 i i i i=1

Td Td is a direct consequence of (45) and properties of η: Define D := i=1 (Ii ×i Vi ) and W := i=1 Wi . Shrinking the sets Ii , Vi , Wi , if necessary, we may assume that f, ∇µ f, η, ∇η, ∆η are uniformly bounded on W × D. Since ψµ = ηf (µ, ·) and ∇µ ψµ = η∇µ f , the integrand in (48) is equal to |∆(ηf )|p |∇µ f |p |η|2−p |f |2−2p . By the product rule, ∆(ηf ) = η ∆f + 2∇η · ∇f + (∆η)f. Hence, for p < 5/4, there exists a constant C < ∞, independent of µ ∈ W , such that on D |∆(ηf )|p |∇µ f |p |η|2−p |f |2−2p ≤ C |∆f |p |∇µ f |p |f |2−2p + |∇f |p |∇µ f |p |f |2−2p + |f |2−p . Pd The first term is controlled by the l = 2 case of (46), since ∆f = j=1 ∂j2 f and ∇µ f is uniformly bounded. The second term is controlled by the l = 1 case of (46), again using the uniform boundedness of ∇µ f . Finally, the third term is uniformly integrable on D, because f is bounded on W × D, D has finite measure, and 2 − p > 0. Therefore Z p ∆ψµ (y) · ∇µ ψµ (y) sup |ψµ (y)|2 dy < ∞, ψµ (y)2 µ∈W D which proves (48). In summary, (48) yields an open covering of Q × K by sets W × V with Z p ∆ψµ (y) · ∇µ ψµ (y) sup |ψµ (y)|2 dy < ∞ ψµ (y)2 µ∈W V By compactness, there is a global constant with Z p ∆ψµ (y) · ∇µ ψµ (y) sup |ψµ (y)|2 dy ≤ C < ∞. 2 ψ (y) µ µ∈Q Ω 5/4

Since H = −∆ + V, this, the compactness of Ω, the analyticity of f and the assumption that V ∈ Lloc imply that there is a global constant with Z p Hψµ (y) · ∇µ ψµ (y) sup |ψµ (y)|2 dy ≤ C < ∞, ψµ (y)2 µ∈Q Ω which yields (35).

PS-Clip-VMC: Provably Convergent Stochastic Optimization

In this section, we propose a robust gradient estimator based on per-sample clipping of both local energies and gradients. We show that, under weak moment assumptions, the resulting optimization method converges to a stationary point both in expectation and with high probability. The estimator and some useful notation are introduced in the following setting. Setting 4.1 (PS-Clip-VMC). Assume Setting 2.1 and consider a parametrized model RP ∋ θ 7→ ψθ ∈ H 2 (Ω). θ ,ψθ ⟩ = E[Eθ ]. For θ ∈ RP define the loss L(θ) := ⟨Hψ ∥ψθ ∥2 L2

Let L∗ := inf θ∈RP L(θ) and assume that L∗ > −∞ (this holds if the bilinear form associated with H is bounded below; see for example [23, Theorem 4.3] for suitable conditions for this to hold). For n ∈ N, α, β ∈ (0, ∞), θ ∈ RP , and i.i.d samples X1 , . . . Xn ∼ Xθ , set Ei = EL,θ (Xi ), Wi = Wθ (Xi ) and let the clipped gradient estimator be defined as n

Gn (θ, α, β) =

2 X 2X γα,p,k (Ek Wk ) − 2 γβ,2,k (Ek ) γβ,2,k (Wj ) , n n −n k=1

(49)

k̸=j

where γα,p,k (·) is the clipping function given by (

αk p γα,p,k (v) := v · min 1, |v|

) for any v ∈ Rd .

(50)

For M, n ∈ N, θ1 ∈ RP , α, β ∈ R, η1 , . . . , ηM −1 ∈ (0, ∞) define the PS-Clip-VMC parameter update rule θm+1 := θm − ηm Gn (θ, α, β),

m = 1, . . . , M − 1,

where for each parameter update θm 7→ θm+1 , the clipped gradient estimator Gn (θ, α, β) is evaluated on i.i.d. samples E1 , . . . , En ∼ Eθm and W1 , . . . , Wn ∼ Wθm . Finally, we let ∆1 := L(θ1 ) − L∗ > 0 the initialization gap. We now state the assumptions under which the convergence results in this section will be proved. Assumptions 4.2. Assume Settings 4.1 and 2.1, and that there are p ∈ (1, 2] and u, σ, C ∈ (0, ∞) such that for any θ ∈ RP , (i) E |Eθ |2 ≤ σ 2 h i 2 (ii) E |Wθ | ≤ σ 2 p

(iii) E [|Eθ Wθ | ] ≤ up (iv) |∇θ L(θ1 ) − ∇θ L(θ2 )| ≤ C|θ1 − θ2 | for any θ1 , θ2 ∈ RP . Remark 4.3. Instead of Item (iv) in Assumptions 4.2, one could also assume h i 2 (v) E |Hθ ψθ /ψθ | ≤ σ 2 h i 2 (vi) E |∇θ Hψθ /ψθ | ≤ σ 2 , which together with (i)-(iii) imply (iv); see [1]. In (v), Hθ ψθ is the Hessian of ψθ . Remark 4.4. It is common practice in variational Monte Carlo methods for solving the Schrödinger equation to clip the local energies during training [20, 19, 8, 10, 14, 12]. Specifically, for some fixed multiplier α ∈ (0, ∞) and samples E1 , . . . , En ∼ Eθt , W1 , . . . , Wn ∼ Wθt , the gradient estimator used for the parameter update θt 7→ θt+1 is given by ! n n 1X 1X G(θt , α) = clip[µ−ασ,µ+ασ] (Ek ) − clip[µ−ασ,µ+ασ] (Ei ) Wk , n n i=1 k=1

Pn Pn where µ = n1 i=1 Ei and σ = n1 i=1 |Ei − µ|. Our analysis suggests that the empirical success of this clipping strategy is rooted in the heavy-tailed nature of the underlying optimization problem. Moreover, our results show that local-energy clipping is sufficient to guarantee convergence in expectation, but not with high probability; see Appendix B. PS-Clip-VMC addresses this issue by additionally clipping the per-sample logarithmic gradients of the wave function. In Section 5, we provide empirical evidence that this method is significantly more robust than the common practice of only clipping the local energy.

4.1

Convergence in Expectation

In this section we prove that under Assumptions 4.2 PS-Clip-VMC converges in expectation to a stationary point. The result is formulated in the following theorem. Theorem 4.5. Assume Settings 4.1 and Assumptions 4.2. Let α = u and β ∈ (0, ∞]. Then, there is a constant Ku,σ,β > 0 that only depends on u, σ and β, such that for any ηm < 1/C the iterates generated by PS-Clip-VMC satisfy M X 2∆1 Ku,σ,β ηm E(|∇θ L(θm )|2 ) ≤ PM + 2(p−1) . PM n p m=1 m=1 ηm m=1 ηm p

In particular, if ηm = 1/(2C) is constant and n ≥ M 2(p−1) we obtain that M

∆1 C + Ku,σ,β 1 X E(|∇θ L(θm )|2 ) ≤ . M m=1 M To prove Theorem 4.5 we begin with the following proposition which provides a standard estimate of the average squared gradient norm of the iterates generated by variational Monte Carlo SGD. Proposition 4.6. Assume Settings 4.1 and Assumptions 4.2(iv). Then, for any ηm < 1/C, the iterates of variational Monte-Carlo SGD satisfy M X

ηm |∇θ L(θm )|2 ≤ 2∆1 +

m=1

M X

ηm |∇θ L(θm ) − Gn (θm , α, β)|2 .

(51)

m=1

Proof. For simplicity we write Gn (θm ) for Gn (θm , α, β). Then, for every m = 1, . . . , M − 1 we have that L(θm+1 ) − L(θm ) ≤ −ηm ⟨∇θ L(θm ), Gn (θm )⟩ +

2 Cηm |Gn (θm )|2 2

2 Cηm ηm |∇θ L(θm )|2 + |Gn (θm )|2 − |∇θ L(θm ) − Gn (θm )|2 + |Gn (θm )|2 2 2 ηm |∇θ L(θm )|2 − |∇θ L(θm ) − Gn (θm )|2 , ≤− 2

=−

where the equality in the second line is simply the polarization identity for inner products. Summing over m and telescoping we obtain that M X m=1

ηm |∇θ L(θm )|2 ≤ 2∆1 +

M X

ηm |∇θ L(θm ) − Gn (θm )|2 .

m=1

By Proposition 4.6, we see that, h to prove convergence of the i algorithm, we need to find an upper bound 2 for the expected squared error E |∇θ L(θm ) − Gn (θm , α, β)| . Note that, since we do not assume |∇θ L(θ)| to have bounded variance, this would not be possible if, instead of Gn , one were to use a non-robust gradient estimator such as the empirical mean. We begin with the following lemma, which provides a generic in-expectation error bound for random variables with a finite p-th moment. In particular, when applied to the random variables Yk = Ek Wk , it yields an estimate of the error generated by the first term in (49). 20

Lemma 4.7. Let Y1 , . . . , Yn be i.i.d. random variables such that E |Y1 | ≤ up < ∞ for some p ∈ (1, 2]. Then, for any α ≥ u we have that n

1X E γα,p,k (Yk ) − E(Y1 ) n

≤ 8α2 n−

2(p−1) p

k=1

where γα,p,k (·) is the clipping function defined in (50). Proof. For simplicity, we set Ȳk = γα,p,k (Yk ) in the following. We set, µ = E(Y1 ) and µk = E(Ȳk ). Then, ! ! 1 αk p |µk − µ| = E − 1 Yk 1(|Yk |>αk1/p ) ≤ E Yk 1(|Yk |>αk1/p ) |Yk | (52) 1−p 1−p 1−p p p 1−p ≤ E |Yk | 1(|Yk |>αk1/p ) α k p ≤ u α k p . Hence, we obtain that n

1X E (µk − µ) n k=1

1 ≤ 2 n

n X

p 1−p

u α

1−p p

≤

k=1

= 4u2p α2(1−p) n

2(p−1) − p

u2p α2(1−p) p1 2 2n n2

≤ 4α2 n

2(p−1) − p

(53)

Further, we have that n

1X Ȳk − µk E n

k=1

≤

1 X 1 X 2 2 E Ȳk − µk ≤ 2 E Ȳk 2 n n 1 n2

k=1 n X

k=1

2−p p

2(p−1) p

≤ 4α2 n−

E |Yk | α2−p k

≤

k=1

p 2−p −

= 4u α

4up 2−p p2 α n n2 2(p−1) p

(54)

where in the second inequality, we used the fact that by definition of γk we have that Ȳk Ȳk

2−p

2−p p

= Ȳk

Ȳk

2−p

≤

. Combining (53) and (54), we obtain that n

1X Ȳk − µ n k=1

≤E

1X Ȳk − µk n

k=1

1X (µk − µ) n

≤ 8α2 n−

2(p−1) p

k=1

which concludes the proof. We now prove a similar error bound for the second term in (49). The proof follows a similar idea to the previous result. However, it requires some additional steps, since, due to the double sum, the random variables involved are not necessarily independent. Lemma 4.8. Assume Settings 4.1 and Assumptions 4.2(i-iii). Then, for any β > 0,   2 ! 2 X 2σ 16σ 2 σ 4 1 3   E γβ,2,k (Ek )γβ,2,k (Wj ) − E(E1 )E(W1 )  ≤ 2 + + 1+ √ . n(n − 1) n−1 β2 n β n k̸=j

Proof. We set Ēk = Ek γβ,2,k (Ek ), W̄k = Wk γβ,2,k (Wk ) and µk,j = E(Ēk W̄j ). First, we note that since Ek and Wj are independent for any k ̸= j we have that     2 2 2 X X  X  E Ēk W̄j − E(E1 )E(W1 )  ≤ E Ēk W̄j − µk,j  + E[Ek Wj ] − µk,j . (55) k̸=j

k̸=j

We begin by bounding the first term on the right-hand side of (55).   2 X  X X E[|Ēk W̄j − µkj |2 ] + E[(Ēk W̄j − µkj )(Ēj W̄k − µjk )] Ēk W̄j − µkj  = E k̸=j

k̸=j

X X j1 ̸=j2 j1 ,j2 ̸=k

E (Ēk W̄j1 − µkj1 )(Ēk W̄j2 − µkj2 )

X X k1 ̸=k2 k1 ,k2 ̸=j

E (Ēk1 W̄j − µk1 j )(Ēk2 W̄j − µk2 j ) E (Ēk1 W̄j1 − µk1 j1 )(Ēk2 W̄j2 − µk2 j2 )

k1 ,k2 ,j1 ,j2 all different

X E(|Ēk W̄k |)E(|Ēj W̄j |) − µkj µjk E(Ēk2 )E(|W̄j |2 ) − |µkj |2 + k̸=j

k̸=j

X X k

≤

(56)

E(Ēk2 )E(W̄j1 )E(W̄j2 ) − µkj1 µkj2

j1 ̸=j2 j1 ,j2 ̸=k

X X j

k1 ̸=k2 k1 ,k2 ̸=j

σ4 +

k̸=j

E(Ēk1 )E(Ēk2 )E(W̄j2 ) − µk1 j µk2 j

2σ 4 +

X X k

σ4 +

X X

j1 ̸=j2 j1 ,j2 ̸=k

σ4

k1 ̸=k2 k1 ,k2 ̸=j

= (3n(n − 1) + 2n(n − 1)(n − 2))σ 4 . Where in the second equality we use that Ēk and W̄j are independent for k ̸= j. In the last inequality in (56) we bounded the terms of the form E(|Ēk W̄k |) by σ 2 . A slightly tighter bound is given by min{u, σ 2 }, but we omit it here to avoid cluttering. For the second term on the right-hand side of (55), we have that X

E(Ek Wj ) − µk,j ≤

k̸=j

X E |Ek Wj |1{|Ek |>β √k , |Wj |≤β √j} k̸=j

X E |Ek Wj |1{|Ek |≤β √k, |Wj |>β √j} k̸=j

X + E |Ek Wj |1{|Ek |>β √k, |Wj |>β √j} k̸=j

≤

σ E |Ek |1{|Ek |>β √k} + σ E |Wj ||Wj |>β √j}

(57)

k̸=j

X E |Ek Wj |{|Ek |>β √k, |Wj |>β √j} k̸=j

X σ E |Ek |2 σ E |Wj |2 E |Ek |2 E |Wj |2 √ √ √ ≤ + + β j β 2 kj β k k̸=j p 3√ √ 4σ n(n − 1) 4σ 3 n(n − 1) 4σ 4 n(n − 1) 2σ ≤ + ≤ 1+ √ , 2 β β β β n where again we use that Ek and Wj are independent for k ̸= j and E(|Ek |), E(|Wj |) < σ. Combining

(55), (56), and (57), we obtain that 

X 1 Ēk W̄j − E(E1 )E(W1 ) n(n − 1)

≤

k̸=j

1 n2 (n − 1)2



 X Ēk W̄j − µk,j  E k̸=j 2

X 1 + 2 E(Ej Wk ) − µk,j 2 n (n − 1) j̸=k

2 (3 + 2(n − 2))n(n − 1)σ 4 2σ 16σ 6 n(n − 1)2 √ + 1 + 2 2 n (n − 1) β 2 n2 (n − 1)2 β n ! 2 3 2σ 16σ 2 σ 4 = 2+ + 1+ √ . n−1 β2 n β n ≤

We can now combine Lemma 4.7 and Lemma 4.8 to obtain an estimate for the expected squared error of the gradient estimator defined in (49). Corollary 4.9. Assume Setting 4.1 and Assumptions 4.2(i-iii). Then, for any β > 0

E |Gn (θ, u, β) − ∇θ L(θ)|

≤ 64u n

2(1−p) p

! 2 2σ 3 16σ 2 8σ 4 + 1+ √ 2+ n−1 β2 n β n

2(1−p) =O n p . Proof. Setting Yk = Ek Wk we have by Lemma 4.7 and Lemma 4.8 that

E |Gn (θ, u, β) − ∇θ L(θ)|

2X (γσ,p,k (Yk ) − E(Y1 )) ≤ 2E n

k=1

X 2 (γβ,2,k (Ek )γβ,2,k (Wj ) − E(E1 )E(W1 )) + 2E n(n − 1) k̸=j ! 2 2 2(1−p) 3 2σ 8σ 4 16σ ≤ 64u2 n p + 2 + + 1+ √ . n−1 β2 n βn

Remark 4.10. Note that if no clipping is applied to the second term in (49), i.e. if “β = ∞”, then the error bound in Corollary 4.9 becomes 4 2(1−p) 3 8σ 2 E |Gn (θ, α, β) − ∇θ L(θ)| ≤ 64u2 n− p + 2 + . n−1 n In particular, convergence in expectation can be established, even without clipping the second term in (49). However, this is no longer the case when considering convergence with high-probability. We can now prove Theorem 4.5 Proof of Theorem 4.5. By Proposition 4.6, we know that the iterates generated by the PS-Clip-VMC satisfy M X m=1

ηm |∇θ L(θm )|2 ≤ 2∆1 +

M X m=1

ηm |∇θ L(θm ) − Gn (θm , α, β)|2 .

(58)

Further, by Corollary 4.9 there is a constant Ku,σ,β > 0, which only depends on u, σ and β, such that E(|∇θ L(θ) − Gn (θ)|2 ≤ Ku,σ,β n

2(1−p) p

), for any θ ∈ Rd . Hence, for any m = 1, . . . , M we have that

2(1−p) E |∇θ L(θm ) − Gn (θm )|2 = E Eθm |∇θ L(θm ) − Gn (θm )|2 ≤ Ku,σ,β n p .

(59)

Plugging (59) into (58), we obtain that M X

ηm E(|∇θ L(θm )|2 ) ≤ 2∆1 + Ku,σ,β n

m=1

and dividing both sides by

4.2

M X

2(1−p) p

ηm ,

m=1

m=1 ηm concludes the proof.

High Probability Convergence

Due to the high computational cost of each training run, it is generally desirable to establish high-probability convergence guarantees for an optimization algorithm in the context of VMC. That is, one seeks to prove that, with probability at least 1−δ, the considered convergence measure depends at most polylogarithmically on log(1/δ). The following theorem provides a high-probability convergence guarantee for PS-Clip-VMC under Assumptions 4.2. Theorem 4.11. Assume Settings 4.1 and Assumptions 4.2. Let δ ∈ (0, 1) and set α = u/(log(4/δ) + 1/4)1/p and β = σ/(log(4/δ) + 1/4)1/2 . Then, for any ηm < 1/C and n ≥ log((4M )/δ) + 1/4, we have that with probability at least 1 − δ the iterates generated by PS-Clip-VMC satisfy 1 PM

M X

m=1 ηm m=1

2∆1

ηm |∇θ L(θm )| ≤ PM

m=1 ηm

2 2

+ 4(7u + 94σ )

log(4M/δ) + 1/4 n

2(p−1) p .

In particular, if ηm = 1/(2C) is constant and n ≥ M 2(p−1) we obtain that M

1 X ∆1 C + 4(7u + 94σ 2 )2 (log(4M/δ) + 1/4) |∇θ L(θm )|2 ≤ M m=1 M

2(p−1) p

Remark 4.12. Note that, in general, n ≫ log((4M )/δ), so that the assumption n ≥ log((4M )/δ) + 1/4 in Theorem 4.11 is not a significant restriction. Similarly as in the previous section, we begin by providing a generic high probability upper bound for random variables with a finite p-th moment. p

Lemma 4.13. Let Y1 , . . . , Yn be i.i.d. random variables such that E |Y1 | ≤ up < ∞ for some p ∈ (1, 2], and α̃ 1/4 for any δ ∈ (0, 1) and α̃ ≥ σ set α = log(1/δ) we have that 1/p . Then, with probability at least 1 − δe n

1X γα,p,k (Yk ) − E(Y1 ) ≤ 7α̃ n k=1

log(1/δ) n

p−1 p .

Proof. For simplicity we set Ȳk = γα̃,p,k (Yk ). The proof of this result follows the same idea as the one of [6, Lemma 1]. First, note that for every k = 1, . . . , n we have that E Ȳk − E(Ȳk ) = 0 , p1 n |Ȳk − E(Ȳk )| ≤ 2α̃ and log(1/δ) 2−p p n 2 p 2−p E |Ȳk − E(Ȳk )| ≤ u α̃ . log(1/δ) 24

Hence, using the second inequality in Lemma A.1 we obtain that with probability at least 1 − δe1/4 s p1 2−p n p X 1 n log(1/δ) n log(1/δ) Ȳk − E(Ȳk ) ≤ 2α̃ + 8up α̃2−p n log(1/δ) n log(1/δ) n k=1 p−1 p−1 log(1/δ) p √ log(1/δ) p = α̃ 2 + 8(uα̃−1 )p/2 < 5α̃ . n n

(60)

Combining this, with (52) we obtain that with probability at least 1 − δe1/4 E(Y1 ) −

n n n 1 X 1X 1 X E(Yk ) − E(Ȳk ) + E(Ȳk ) − Yk γk Ȳk ≤ n n n k=1 k=1 k=1 p−1 n p p−1 1−p 1 X p 1−p log(1/δ) ≤ u α̃ log(1/δ) p k p + 5α̃ n n k=1 p−1 p p−1 1 2 p 1−p log(1/δ) log(1/δ) p n p + 5α̃ ≤ u α̃ n n p−1 p p−1 1−p log(1/δ) p 1−p ≤ 2u α̃ log(1/δ) p n p + 5α̃ n p−1 log(1/δ) p ≤ 7α̃ , n

where in the second inequality we used (52) to estimate the first term and (60) for the second one. We can now provide a high probability error bound for the gradient estimator defined in Setting 4.1. Lemma 4.14. Assume Settings 4.1 and Assumptions 4.2(i-iii). For any δ ∈ (0, 1), set α = u/ log(4/δ)1/p and β = σ/ log(4/δ)1/2 . Then, with probability at least 1 − e1/4 δ r p−1 log(4/δ) p log(4/δ) log(4/δ) 12σ 2 2 2 + 24σ + 58σ + n n n n ! p−1 log(4/δ) p =O . n

1 |Gn (θ, α, β) − ∇θ L(θ)| ≤ 7u 2

Proof. Let Yk = Ek Wk and define the clipped random variables Ēk = γβ,2,k (Ek ), W̄k = γβ,2,k (Wk ) and Ȳk = γα,p,k (Yk ). Then, with this notation we have that  !  n X X 1 1 1 (Gn (θ, α, β) − ∇θ L(θ)) = Ȳk − E(Y1 ) + E(E1 )E(W1 ) − 2 Ēk W̄j  . (61) 2 n n −n k=1

k̸=j

Let ρ = δ/4, so that α = u/ log(1/ρ)1/p and β = σ/ log(1/ρ)1/2 . Then, for the first term on the right-hand side of (61), we know by Lemma 4.13 that with probability at least 1 − e1/4 ρ n

1X Yk γα,p,k (Yk ) − E(Y ) ≤ 7u n k=1

log(1/ρ) n

p−1 p .

(62)

To estimate the second term in (61) we rewrite it as follows: E[E1 ]E[W1 ] −

X X 1 1 Ēj W̄k = E[E1 ]E[W1 ] − E[Ēj ]E[W̄k ] n(n − 1) n(n − 1) j̸=k j̸=k | {z } I

X 1 − Ēj W̄k − E[Ēj ]E[W̄k ] n(n − 1) j,k {z } |

(63)

X 1 + Ēj W̄j − E[Z̄j ]E[W̄j ] . n(n − 1) j | {z } III

Note that by independence, we have that E[E1 ]E[W1 ] = E(Ej Wk ) for any j ̸= k. Hence, using (57) we obtain the following bound for the term I on the right-hand side of (63):

X 1 E[Ej Wk ] − E[Ēj W̄k ] ≤ (I) = n(n − 1) j̸=k

2σ 1+ √ β n

4σ 3 √ = β n

1+2

log(1/ρ) n

n X

log(1/ρ) 2 4σ . (64) n

To estimate (II) we rewrite the sum as n X

Ēj W̄k − E[Ēj W̄k ] =

j,k=1

n X

E[W̄k ]

n X

(Ēj − E[Ēj ]) +

(Ēj − E[Ēj ])

E[Ēj ]

j=1

k=1 n X j=1

n X

(W̄k − E[W̄k ])

k=1

(W̄k − E[W̄k ]) .

k=1

Note that the random variables Ēj − E[Ēj ] and W̄j − E[W̄j ] satisfy

E Ēj − E(Ēj ) = 0 ,

E W̄j − E(W̄j ) = 0 ,

|Ēj − E(Ēj )| ≤ 2σ

|W̄j − E(W̄j )| ≤ 2σ

n log(1/ρ)

E |Ēj − E(Ēj )|2 ≤ σ 2

E |W̄j − E(W̄j )|2 ≤ σ 2

n log(1/ρ)

Hence, using the same steps as in (60), we obtain that with probability 1 − 2e1/4 ρ 1 (II) = n(n − 1)

n X k=1

E(W̄k )

n X

(Ēj − E(Ēj )) +

j=1

n n X X + (Ēj − E(Ēj )) (W̄k − E(W̄k )) j=1

n X

E(Ēj )

n X

(W̄k − E(W̄k ))

k=1



 12 X 12 n n X 1 log(1/ρ) log(1/ρ)   ≤ E(|W̄k |)5σ + E(|Ēj |)5σ (n − 1) n n j=1 k=1 n2 log(1/ρ) 2 25σ + n(n − 1) n 1/2 log(1/ρ) log(1/ρ) ≤ 20σ 2 + 50σ 2 . n n

(65)

To estimate (III), we note that the term can be rewritten as follows (III) =

n n X X 1 1 Ēj W̄j − E(Ēj W̄j ) + E(Ēj W̄j ) − E(Ēj )E(W̄j ). n(n − 1) j=1 n(n − 1) j=1

(66)

The second term in (67) is deterministic and easily seen to be bounded by σ 2 /(n − 1) ≤ 2σ 2 /n. For the first term, we note that it satisfies E Ēj W̄j − E(Ēj W̄j ) = 0 n 2 |Ēj W̄j − E(Ēj W̄j )| ≤ 2σ log(1/ρ) n 2 4 E |Ēj W̄j − E(Ēj W̄j )| ≤ σ log(1/ρ) Hence, similarly as before, we obtain that with probability at least 1 − ρe1/4 n X 1 2 √ 2 10σ 2 1 Ēj W̄j − E(Ēj W̄j ) ≤ 2σ + 8σ ≤ n(n − 1) j=1 n−1 n

(67)

Hence, plugging (62), (64), (65), and (67) into (61) we obtain that with probability 1 − 4e1/4 ρ !r r p−1 log(1/ρ) log(1/ρ) 2 1 log(1/ρ) p + 1+2 (Gn (θ, α, β) − ∇θ L(θ)) ≤ 7u 4σ 2 n n n r log(1/ρ) log(1/ρ) 12σ 2 2 2 + 50σ + 20σ + n n n r p−1 log(1/ρ) p log(1/ρ) 12σ 2 log(1/ρ) 2 2 = 7u + 24σ + 58σ + n n n n The statement in the Lemma follows by recalling that ρ = δ/4. Corollary 4.15. Let θ1 , . . . , θM be arbitrary weights, δ ∈ (0, 1) and n ≥ log((4M )/δ) + 1/4. Set α = u/(log(4/δ) + 1/4)1/p and β = σ/(log(4/δ) + 1/4)1/2 . Then, under the same conditions of Lemma 4.14, we have that p−1 log(4M/δ) + 1/4 p max |Gn (θm , α, β) − ∇θ L(θ)| ≤ 2(7u + 94σ 2 ) . m=1,...,M n Proof. Follows immediately by applying union bound to Lemma 4.14. The proof of Theorem 4.11 now follows using similar steps to those of its in Expectation counterpart. Proof of Theorem 4.11. By Proposition 4.6 we know that the iterates of the algorithm satisfy M X

ηm |∇θ L(θm )|2 ≤ 2∆1 +

m=1

M X

ηm |∇θ L(θm ) − Gn (θm , α, β)|2 .

m=1

Hence, applying Corollary 4.15 we obtain that with probability at least 1 − δ M X

2 2

ηm |∇θ L(θm )| ≤ 2∆1 + 4(7u + 94σ )

m=1

and the result follows by dividing both sides by

m=1 ηm .

log(4M/δ) + 1/4 n

2(p−1) M p X m=1

ηm ,

Experiments

We empirically tested the robustness of PS-Clip-VMC and compared its performance to the standard practice of only clipping the local energy. To this end, we use the FermiNet implementation from [15] with 256 hidden dimensions and 16 Slater determinants. Although our theoretical results are based on a clipping threshold which increases with each sample, we found that a constant clipping threshold, analogue to the one commonly used for local energy clipping was more effective. Specifically, let n ∈ N be the batch size (in our case 2048), E1θ , . . . , Enθ be samples ofP the local n 1 θ θ θ energy and W1θ , . . . , Wnθ corresponding samples of ∇ log |ψ |. Further, let µ = µ(E , . . . , E ) = θ θ n 1 i=1 Ei , n Pn Pn 1 1 θ ′ ′ θ θ ′ ′ θ θ θ θ θ σP = σ(E1 , . . . , En ) = n i=1 |Ei − µ| and µ = µ (W1 , . . . , Wn ) = n i=1 Ei , σ = σ (W1 , . . . , Wn ) = n 1 θ ′ i=1 |Wi − µ |. Then, the energy gradient estimator for PS-Clip-VMC used in our experiments is given n by ) ( ! n n 1X 1X µ′ + βσ ′ θ θ Wiθ , (68) G(θ, α, β) = min 1, clip[µ−ασ,µ+ασ] Ei clip[µ−ασ,µ+ασ] Ei − n i=1 n i=1 Wiθ where

  a if x < a clip[a,b] (x) = x if a ≤ x ≤ b .   b if x > b

For both methods, we use the default hyperparameters from the FermiNet GitHub repository, except for the number of decorrelation MCMC steps, which we increase to 30. The model was trained for 2 · 105 iterations with a batch size of 2048; note that this corresponds to half of the batch size normally used to train FermiNet. We refer to Table 2 for a detailed list of the hyperparameters used. To efficiently compute per-sample gradient norms, we use Google’s JAX Privacy API (https://github.com/google-deepmind/jax_privacy). Table 1 compares the final energies obtained using only standard (centered) local energy clipping with those obtained using the additional PS-gradient clipping. We observe that PS-Clip-VMC outperforms the standard method in both instances, with the performance difference for Argon being particularly significant. The poor performance of the standard method on Argon can be explained by examining Figure 1. While PSClip-VMC converges steadily, the training trajectory of the standard method is more unstable. In particular, it exhibits a cusp around step 55,000, after which it is unable to recover a lower energy. Parameter

Value

Training

Learning rate at time t Initial learning rate Learning rate decay Local energy clipping PS-Gradient clipping

lr0 (1 + t/t0 )−1 lr0 = 0.05 t0 = 105 α = 5.0 β = 5.0 or “β = ∞”

Pretraining

Pretraining optimizer Pretraining iterations Pretraining basis set

Adam 104 ccpvdz

MCMC

Batch size Decorrelation steps

2048 30

KFAC

Norm constraint Damping

10−3 10−3

Table 2: Default hyperparameters used to train the network. β = ∞ corresponds to the standard practice of only clipping the local energy. See also (68) for the definition of the parameters α and β.

Conclusion

In this work, we demonstrated that the stochastic optimization problem arising in variational Monte Carlo for the electronic Schrödinger equation is naturally heavy-tailed. Specifically, we proved that for commonly used ansatz classes, such as Slater–Jastrow wave functions with variable-exponent Slater-type orbitals, the random variable used to estimate the stochastic gradient generally does not possess moments of order greater than 3/2. On the other hand, we showed that for any analytic, compactly supported ansatz, this random variable possesses moments of order at least 5/4. Motivated by these observations, we propose a robust optimization algorithm that guarantees convergence of VMC both in expectation and with high probability in the weak-moment regime. We validate our theoretical results through numerical experiments using FermiNet.

References [1]

Nilin Abrahamsen et al. “Convergence of variational Monte Carlo simulation and scale-invariant pre-training”. In: Journal of Computational Physics 513 (2024), p. 113140. doi: 10 . 1016 / j . jcp . 2024.113140.

[2]

Shmuel Agmon. “Bounds on exponential decay of eigenfunctions of Schrödinger operators”. In: Schrödinger Operators: Lectures given at the 2nd 1984 Session of the Centro Internationale Matematico Estivo (CIME) held at Como, Italy, Aug. 26–Sept. 4, 1984. Springer, 2006, pp. 1–38.

[3]

JM Aldaz, Ognyan Kounchev, and Hermann Render. “Bernstein operators for exponential polynomials”. In: Constructive Approximation 29.3 (2009), pp. 345–367.

[4]

Victor Armegioiu et al. “Functional Neural Wavefunction Optimization”. In: arXiv preprint arXiv:2507.10835 (2025).

[5]

Federico Becca and Sandro Sorella. Quantum Monte Carlo approaches for correlated systems. Cambridge University Press, 2017.

[6]

Sebastien Bubeck, Nicolo Cesa-Bianchi, and Gabor Lugosi. “Bandits With Heavy Tail”. In: IEEE Trans. Inf. Theor. 59.11 (2013), pp. 7711–7717.

[7]

Subhas J. Chakravorty et al. “Ground-state correlation energies for atomic ions with 3 to 18 electrons”. In: Physical Review A 47.5 (1993), pp. 3649–3670.

[8]

Nicholas Gao and Stephan Günnemann. “Ab-Initio Potential Energy Surfaces by Pairing GNNs with Neural Wave Functions”. In: International Conference on Learning Representations. 2022. doi: 10. 48550/arXiv.2110.05064. arXiv: 2110.05064.

[9]

Nicholas Gao and Stephan Günnemann. “Generalizing Neural Wave Functions”. In: Proceedings of the 40th International Conference on Machine Learning. Vol. 202. Proceedings of Machine Learning Research. PMLR, 2023, pp. 10708–10726. doi: 10.48550/arXiv.2302.04168. arXiv: 2302.04168.

[10]

Nicholas Gao and Stephan Günnemann. “Sampling-free Inference for Ab-Initio Potential Energy Surface Networks”. In: The Eleventh International Conference on Learning Representations. 2023. doi: 10.48550/arXiv.2205.14962. arXiv: 2205.14962.

[11]

Leon Gerard, Philipp Grohs, and Michael Scherbela. “Deep learning variational Monte Carlo for solving the electronic Schrödinger equation”. In: Numerical Analysis Meets Machine Learning. Vol. 25. Handbook of Numerical Analysis. North-Holland, 2024, pp. 231–292. doi: 10.1016/bs.hna.2024.05.010.

[12]

Leon Gerard et al. “Gold-standard solutions to the Schrödinger equation using deep learning: How much physics do we need?” In: Advances in Neural Information Processing Systems. Vol. 35. 2022, pp. 10282–10294. doi: 10.48550/arXiv.2205.09438. arXiv: 2205.09438.

[13]

Leon Gerard et al. “Transferable neural wavefunctions for solids”. In: Nature Computational Science 5 (2025), pp. 1147–1157. doi: 10.1038/s43588-025-00872-z. arXiv: 2405.07599.

[14]

Ingrid von Glehn, James S. Spencer, and David Pfau. “A Self-Attention Ansatz for Ab-initio Quantum Chemistry”. In: The Eleventh International Conference on Learning Representations. 2023. doi: 10. 48550/arXiv.2211.13672. arXiv: 2211.13672.

[15]

Google DeepMind. FermiNet. GitHub repository. 2020. url: https://github.com/google-deepmind/ ferminet.

[16]

Philipp Grohs and Martin Rathmair. “Stable Gabor phase retrieval for multivariate functions”. In: Journal of the European Mathematical Society 24.5 (2021), pp. 1593–1615.

[17]

Stephen J Gustafson et al. Mathematical concepts of quantum mechanics. Vol. 33. Springer, 2003.

[18]

Trygve Helgaker, Poul Jorgensen, and Jeppe Olsen. Molecular electronic-structure theory. John Wiley & Sons, 2013.

[19]

Jan Hermann, Zeno Schätzle, and Frank Noé. “Deep-neural-network solution of the electronic Schrödinger equation”. In: Nature Chemistry 12.10 (2020), pp. 891–897.

[20]

Jan Hermann et al. “Ab initio quantum chemistry with neural-network wavefunctions”. In: Nature Reviews Chemistry 7 (2023), pp. 692–709. doi: 10.1038/s41570-023-00516-8. 30

[21]

Morris W Hirsch. Differential topology. Springer Science & Business Media, 2012.

[22]

Tosio Kato. “On the eigenfunctions of many-particle systems in quantum mechanics”. In: Communications on Pure and Applied Mathematics 10.2 (1957), pp. 151–177.

[23]

Tosio Kato. Perturbation theory for linear operators. Second. Grundlehren der Mathematischen Wissenschaften, Band 132. Springer-Verlag, Berlin-New York, 1976, pp. xxi+619.

[24]

Jonas Moritz Kohler and Aurelien Lucchi. “Sub-sampled cubic regularization for non-convex optimization”. In: Proceedings of the 34th International Conference on Machine Learning - Volume 70. ICML’17. 2017, pp. 1895–1904.

[25]

Jacob Korevaar and Jan Wiegerinck. Several complex variables. Korteweg-de Vries Institute for Mathematics Amsterdam, 2017.

[26]

Tianyou Li et al. Convergence Analysis of Stochastic Gradient Descent with MCMC Estimators. 2024. arXiv: 2303.10599.

[27]

Jeffmin Lin, Gil Goldshlager, and Lin Lin. “Explicitly antisymmetrized neural network layers for variational Monte Carlo simulation”. In: Journal of Computational Physics 474 (2023), p. 111765.

[28]

Bernard Malgrange. “The preparation theorem for differentiable functions”. In: Differential Analysis, Bombay Colloq. 1964, pp. 203–208.

[29]

Davide Nobile and Philipp Grohs. “Robust and Fast Training via Per-Sample Clipping”. In: arXiv preprint arXiv:2605.02701 (2026).

[30]

Tianyu Pang, Shuicheng Yan, and Min Lin. O(N 2 ) Universal Antisymmetry in Fermionic Neural Networks. 2022. doi: 10.48550/arXiv.2205.13205. arXiv: 2205.13205.

[31]

David Pfau et al. “Ab initio solution of the many-electron Schrödinger equation with deep neural networks”. In: Physical review research 2.3 (2020), p. 033429.

[32]

Allan Pinkus. N-widths in Approximation Theory. Springer Science & Business Media, 2012.

[33]

R Michael Range. Holomorphic functions and integral representations in several complex variables. Vol. 108. Springer Science & Business Media, 1998.

[34]

Nicolas Renaud. “QMCTorch: Molecular Wave Function with Neural Components for Energy and Force Calculations”. In: Methods 4 (2025), p. 4.

[35]

Michael Scherbela, Leon Gerard, and Philipp Grohs. Towards a Foundation Model for Neural Network Wavefunctions. Preprint title/version of the transferable molecular neural-wavefunction work. 2023. doi: 10.48550/arXiv.2303.09949. arXiv: 2303.09949 [physics.chem-ph].

[36]

Michael Scherbela, Leon Gerard, and Philipp Grohs. “Towards a transferable fermionic neural wavefunction for molecules”. In: Nature Communications 15 (2024), p. 120. doi: 10.1038/s41467-02344216-9.

[37]

Michael Scherbela, Leon Gerard, and Philipp Grohs. “Variational Monte Carlo on a Budget: Finetuning pre-trained Neural Wavefunctions”. In: Advances in Neural Information Processing Systems. Vol. 36. 2023. doi: 10.48550/arXiv.2307.09337. arXiv: 2307.09337.

[38]

Michael Scherbela et al. Accurate Ab-initio Neural-network Solutions to Large-Scale Electronic Structure Problems. 2025. arXiv: 2504.06087 [physics.comp-ph].

[39]

Michael Scherbela et al. “Solving the electronic Schrödinger equation for multiple nuclear geometries with weight-sharing deep neural networks”. In: Nature Computational Science 2.5 (2022), pp. 331–341. doi: 10.1038/s43588-022-00228-x. arXiv: 2105.08351.

[40]

James S. Spencer et al. Better, Faster Fermionic Neural Networks. 2020. doi: 10.48550/arXiv.2011. 07125. arXiv: 2011.07125 [physics.chem-ph].

[41]

Attila Szabo and Neil S Ostlund. Modern quantum chemistry: introduction to advanced electronic structure theory. Courier Corporation, 2012.

[42]

Gerald Teschl. Mathematical methods in quantum mechanics. Vol. 157. American Mathematical Soc., 2014.

[43]

Julien Toulouse, Roland Assaraf, and Cyrus J Umrigar. “Introduction to the variational and diffusion Monte Carlo methods”. In: Advances in quantum chemistry. Vol. 73. Elsevier, 2016, pp. 285–314.

[44]

JR Trail. “Heavy-tailed random error in quantum Monte Carlo”. In: Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 77.1 (2008), p. 016703.

[45]

Holger Wendland. Scattered data approximation. Vol. 17. Cambridge university press, 2004.

Additional Technical Results

Lemma A.1 (Vector Bernstein Inequality, see [24, Lemma 18]). Let X1 , . . . , Xn be independent vector-valued random variables with common dimension d, satisfying E[Xi ] = 0 , |Xi | ≤ c and E |Xi |2 ≤ σ 2 , for all i = 1, . . . , n and some c, σ > 0. Then, for every ϵ > 0 we have that ! n 1X −nϵ2 1 P Xi ≥ ϵ ≤ exp + . n i=1 8σ 2 + cϵ 4 In particular, this implies that for any δ ∈ (0, 1) n

c log(1/δ) 1X Xi ≥ + n i=1 n

8σ 2 log(1/δ) n

! ≤ δe1/4 .

In Expectation Convergence of Local Energy Clipping

As mentioned in Remark 4.4, it is common practice in variational Monte Carlo to clip the local energy during training. In this section, we show, using arguments similar to those in Section 4.1, that clipping only the local energy is sufficient to guarantee convergence in expectation. However, unlike PS-Clip-VMC, this approach does not guarantee convergence with high probability. Setting B.1 (VMC with Local Energy Clipping). Assume Setting 2.1 and consider a parametrized model θ ,ψθ ⟩ RP ∋ θ 7→ ψθ ∈ H 2 (Ω). For θ ∈ RP define the loss L(θ) := ⟨Hψ = E[Eθ ]. ∥ψθ ∥2 L2

Let L∗ := inf θ∈RP L(θ) and assume that L∗ > −∞. For n ∈ N, α ∈ (0, ∞), θ ∈ RP , and i.i.d samples X1 , . . . Xn ∼ Xθ , set Ei = EL,θ (Xi ), Wi = Wθ (Xi ) and let the clipped energy gradient estimator be defined as   n n 1 X 1X Gn (θ) = γu,4−p,k (Ek ) − γu,4−p,j (Ej ) Wk , (69) n−1 n j=1 k=1

where γα,p,k (·) is the clipping function defined in (50). For M, n ∈ N, θ1 ∈ RP , α, β ∈ R, η1 , . . . , ηM −1 ∈ (0, ∞) define the Energy-Clip-VMC parameter update rule θm+1 := θm − ηm Gn (θ, α),

m = 1, . . . , M − 1,

where for each parameter update θm 7→ θm+1 , the clipped gradient estimator Gn (θ, α) is evaluated on i.i.d. samples E1 , . . . , En ∼ Eθm and W1 , . . . , Wn ∼ Wθm . As before, ∆1 := L(θ1 ) − L∗ > 0 denotes the initialization gap. Assumptions B.2. Assume Settings 2.1 and B.1, and that there are p ∈ (1, 2] and u, σ, C ∈ (0, ∞) such that for any θ ∈ RP , (i) E |Eθ |2 ≤ σ 2 h i 2 (ii) E |Wθ | ≤ σ 2 32

h i h i 2−p 2 3−p (iii) E |Eθ | · |Wθ | , E |Eθ | |Wθ | ≤ u4−p (iv) |∇θ L(θ1 ) − ∇θ L(θ2 )| ≤ C|θ1 − θ2 | for any θ1 , θ2 ∈ RP . Remark B.3. Using the same arguments as in Section 3, it is straightforward to show that the class of functions considered in Theorem 3.2 also satisfy Assumptions B.2(iii) for any p ∈ ( 32 , 2]. Remark B.4. While, as we will see in the following, local energy clipping guarantees convergence in expectation, it does not provide a high-probability error bound analogous to the one in Lemma 4.14. This can be seen from the fact that, for a fixed θ ∈ RP , any bound of the form P (|Gn (θ) − ∇θ L(θ)| > δ) ≤ ae−bδ

in particular implies that E [|Gn (θ)|p ] < ∞ for any p > 0. However, by the same arguments used in Section 2, Gn (θ) does not, in general, even possess a third moment. Theorem B.5. Assume Settings B.1 and Assumptions B.2. Then, there is a constant Ku,σ > 0 that only depends on u and σ, such that any ηm < 1/C M X 2∆1 Ku,σ ηm E(|∇θ L(θm )|2 ) ≤ PM + 2(2−p) . PM n 4−p m=1 m=1 ηm m=1 ηm 4−p

In particular, if ηm = 1/(2C) is constant and n = M 2(2−p) we obtain that M

1 X ∆1 C + Ku,σ E(|∇θ L(θm )|2 ) ≤ ;. M m=1 M The proof of Theorem B.5 follows the same steps as the one of Theorem 4.5. We begin with the following proposition. Proposition B.6. Assume Settings B.1 and Assumption B.2(iii). Then,   2 n X 2(p−2) 1 E γu,4−p,k (Ek )Wk − E(E1 W1 )  ≤ 8u4 n 4−p . n k=1

Proof. For ease of notation, we set Ēk = Ek γu,4−p,k (Ek ) and define µ = E[E1 W1 ] and µk = E(Ēk Wk ). Then, we have that p−2 p−2 |µ − µk | = |E Ek − Ēk Wk | ≤ E(|Ek 1(|Ek |>u√k) Wk |) ≤ E |Ek |3−p |Wk | up−2 k 4−p ≤ u2 k 4−p . Hence, we obtain that     2 2 n n 2 X X 2(p−2) p−2 2 1 1 1 E (µk − µ)  ≤ 2 E  u2 k 4−p  ≤ 2 2u2 n 4−p = 4u4 n 4−p . n n n k=1

Further, we have that   2 n n n X p 1 1 X 1 X 2   E Ēk − µk ≤ 2 E |Ēk Wk | ≤ 2 E |Ek |2−p |Wk |2 up k 4−p n n n k=1

(70)

k=1

(71)

2(p−2) 4 1 ≤ 2 4u4 n 4−p ≤ 4u4 n 4−p . n

Combining (70) and (71) we obtain that       2 2 2 n n n X X X 2(p−2) 1 1 1 E Ēk Wk − E(E1 W1 )  ≤ E  (µk − µ)  + E  Ēk − µk  ≤ 8u4 n 4−p . n n n k=1

k=1

(72)

Next, we provide an error bound for the double sum term in (??). Proposition B.7. Assume Settings B.1 and Assumptions B.2(i-iii). Then,   2 4 n X 2(p−2) 1 3 σ   E γu,4−p,k (Ek )Wj − E(E1 )E(W1 )  ≤ +2 + 4u4 σ 2 n 4−p . n(n − 1) n−1 n k̸=j

Proof. We again set Ēk = γu,4−p,k (Ek ) and define µ = E[E1 W1 ] and µk = E(Ēk Wk ). Then, we have that     2 2   X X Ēk Wj − E(Ēk )E(Wj )  Ēk Wj − E(E1 )E(W1 )  ≤ E  E k̸=j

k̸=j

 2 X   (E(Ēk ) − E(Ek ))E(Wj )  . + E 

(73)

k̸=j

The first term on the right-hand side of (73) can be estimated using the exact same steps as in (56) to obtain that   2

 X E Ēk Wj − E(Ēk )E(Wj )  ≤ (3n(n − 1) + 2n(n − 1)(n − 2))σ 4 .

(74)

k̸=j

Further, for the second term in (73) we have that     2 2 n X X 4   E(Ēk ) − E(Ek ))  ≤ (n − 1)2 σ 2 4u4 n 4−p , (E(Ēk ) − E(Ek ))E(Wj )  ≤ E  nσ E

(75)

k=1

k̸=j

where the last inequality follows by (71). Combining (74) and (75) we obtain that   2 4 X 1 σ 2 4u4 n 4−p  (3 + 2(n − 2))n(n − 1)σ 4  E Ēk Wj − E(E1 )E(W1 )  ≤ + n(n − 1) n2 (n − 1)2 n2 k̸=j

≤

3 +2 n−1

2(p−2) σ4 + 4u4 σ 2 n 4−p n

Combining Proposition B.6 and Proposition B.7 it is now straightforward to prove the following Lemma. Lemma B.8. Assume Settings B.1 and Assumptions B.2(i)-(iii). Then,     2 n n 2(p−2) X X  1  γu,4−p,k (Ek ) − 1 E γu,4−p,j (Ej ) Wk − E (E1 − E(E1 ))W1  = O n 4−p . n−1 n j=1 k=1

Proof. Setting Ēk = γu,4−p,k (Ek ), we obtain by Propositions B.6 and B.7 that       2 2 n n n X X X 1 1 1   Ēk − E Ēj  Wk − E (E1 − E(E1 ))W1  ≤ 2E  Ēk Wk − E(E1 W1 )  n−1 n j=1 n k=1

k=1



 2 (76) X   + 2E  Ēk Wj − E(E1 )E(W1 )  k̸=j

=O n 34

2(p−2) 4−p

The proof of Theorem B.5 now follows by performing the same steps as in the proof of Theorem 4.5 with the error bound given by Lemma B.8.

Record · ID 307008 · SHA-256 c1bc025dfd087ead
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.