Stat inefficiency

[Pablo_Lopez_Rios]

The sweep does look like a quadric family: over the frames the surface goes from a closed ellipsoid, opens up, passes through something that looks like a hyperboloid, and closes back into an ellipsoid. That is what a family of second-order surfaces does, so at least in this region the node looks like a quadric in the scanned electron's coordinates rather than a curved plane.

Indeed, and there is probably a simple reason for the shapes - the trial wave function used in the plot is at heart a short-ish multideterminant expansion after all. What I tried to convey with the animation though is that a smooth, continuous backflow coordinate transformation x_i(R)=r_i+xi_i(R) cannot possibly take the Hartree-Fock node (a sphere) and make it go through this type of deformation because this would involve mapping points at infinity onto the surface of the sphere at some point in the animation, implying that xi_i(R) should diverge [i.e., so that Psi_HF is evaluated on the sphere (|xi_i|=r_sphere) when r_i is on the opened-up nodal surface at a point infinitely far from the origin].

[Mike Towler]

Vladimir's mechanism is a prediction, and it seems I've been accidentally
testing it. Neil and Pablo have a draft of my new DMC paper since last week;
the sequel derives, on solvable models, an exact propagator for a walker near a
node and an exact, bounded branching factor to go with it, and the sequel to
that puts them into CASINO. They are in my private version, verified against
exact transfer operators at the 1e-5 level, and a registered beryllium ladder
on the Hartree-Fock node has been running on my crappy computer since the
weekend in two arms that differ only in the near-node branching: one with the
usual limdmc-bounded exponent, one with the exact bounded ratio. If the
population-correlation warning comes from near-node excursions, the second
arm's spike fraction and population-weight autocorrelation fall; if it comes
from elsewhere, they don't. We'll see in a day or two.

Meanwhile, one thing about the code. On a determinant-only beryllium trial
under limdmc 4 the population lives about 3 per cent above target, all run,
both arms, and the growth energy is then exactly the mixed energy minus
ln(W/W_target): a limiter-driven population bias that scales as dt to about
0.63, and that a dt times N_w = const sequence cancels rather than removes. And
for calibration: on helium 2 3S, where the node is exact, the default CASINO
step is exact to 1e-5 across dt 0.001 to 0.02, which the DMC paper predicts.
Where the node is right, the code is already right; the question, as Dario's
figure says, is what happens on the sheets and at the crossing, and that is
what the instruments are for.

The follow-up puts limdmc 5 and 6 beside 4 at one rung, since on a bare
determinant the limiter, not the node, owns the linear term; if one of them is
what you would call production now, say so and it goes in the registration.

[Vladimir_Konjkov]

Mike — limdmc 4 here: the CASINO default, unset in every one of my 1673 local DMC directories, and the only scheme PyCasino implements. Put 4 in the registration.

What it is for. I want to optimize the backflow of a single-determinant Be-atom not by emin but on the nodal surface integral of Mitas and Annaberdiyev,
arXiv:2109.01734 — their Eq. (18), the fixed-node energy written as an integral over the node itself rather than as ⟨H⟩. Then compare the fixed-node energies of the two backflows, same determinant, same topology, only the functional they were optimized against differing. My reference point is the emin backflow on the HF node, which sits 2.379 ± 0.054 mHa above the exact non-relativistic energy, -14.667356.

So the comparison has to resolve a few tenths of a mHa between two wave functions of one system, and I would rather get the DMC protocol right the first time. Is dt·N_w = const at limdmc 4 still what you would recommend for that, or does the near-node propagator change what the best sequence is?

[Vladimir_Konjkov]

Bosonic chemistry, anyone?

A small confession from the nodal-surface corner. Eq. (20) of Mitas & Annaberdiyev is beautiful: the fixed-node energy is the bosonic ground-state energy plus one surface integral over the node, weighted by the exact bosonic ground state Φ_B. No curvature, no projection, just an integral.

Then you ask where Φ_B comes from. It is nodeless, so DMC is sign-problem-free for it — but DMC gives you an energy and a mixed distribution, and a Metropolis weight wants the function pointwise. So you write it down the only way anyone ever writes a wave function down: a product of one orbital times a Jastrow, and you optimize it. At which point you notice you have just built a machine for computing the ground-state energy of beryllium as if electrons were bosons.

Which, for any two-electron ground state, is the same number you started with. Congratulations: He and H₂ come out exactly right in the parallel universe too.

I am not complaining — the identity is exact, and the affordable version of §4.1 does earn its keep. But it does mean that anyone optimizing a node this way has to produce a bosonic ground state along the way, at least approximately, and that by-product strikes me as worth a look in its own right. Not for practical reasons: there is no shell structure, no valence and no periodicity to speak of, so the resulting table has one row and the discussion section writes itself. Out of curiosity, then. Sign-problem-free benchmarks for a whole row of the periodic table are cheap to produce and nobody seems to have bothered. Has anyone here?

By the way, in the Jastrow factor for bosons, the electron-electron cusp condition is the same as for antiparallel spins (opposite spins), not parallel spins.

If electrons were bosons, their total wavefunction would have to be completely symmetric under particle exchange. Without a spin degree of freedom to provide antisymmetry, the spatial part of the wavefunction must be strictly symmetric.

In the fermionic world, a symmetric spatial wavefunction corresponds exactly to the singlet state (antiparallel spins). Therefore, the behavior of the bosonic wavefunction at the electron-electron coalescence point (r→0) is mathematically identical to that of antiparallel fermions.

[Mike Towler]

Hi Vladimir, two answers and one result.

Protocol: for a few tenths of a mHa between two wavefunctions of one system, not dt·N_w = const. The constant that cancels the time-step term against the population term belongs to each wavefunction, so what cancels for one backflow does not cancel for the other. Fix the population, run the same dt ladder for both (four rungs), fit each as E0 + c1√dt + c2·dt and look at c1: on a Jastrow-free beryllium determinant I have just found a √dt term of 0.075 Ha au^-1/2 at eighteen standard errors, which no linear trick survives. With a Jastrow it is probably gone, but check. Then plot the difference against dt; if it is flat, its extrapolation is safe whatever the individual curvatures do.

The result I promised: the bounded near-node branching changed nothing in the population statistics on the Hartree-Fock node of beryllium. Your spike fraction was 5.2 to 5.9e-4 with and without it at five time steps from 0.00025 to 0.005, the measured inefficiency 0.95 to 1.16, no difference in either direction (Pablo's caveat on the measured number stands). On a bare determinant the spikes are not near-node excursions; the coalescence tail of a cuspless trial is the candidate, and your Jastrow-and-backflow runs may be a different regime. There is a reason, and it bears on your project: for a determinant of radial orbitals the node is the exchange surface r_i = r_j, and the Laplacian of an antisymmetric function of two radii vanishes there with the function, so the local energy has no 1/s pole on the node, however wrong its topology. The local residual is zero where your Mitas-Annaberdiyev integral is ten mHa; they measure different things. Backflow moves the node off the exchange surface and gives it a pole, which is where your trials and my propagator meet.

Bosonic chemistry: not that I know of, and it is cheap in CASINO (a nodeless product-times-Jastrow trial, DMC exact up to the time step). Your two-electron remark is the whole trap in one line. If you run the row, I would read it.

[Vladimir_Konjkov]

Plan: does Eq. (18) find the DMC optimal node?

Be, CASSCF(2,4)/ano-pVDZ, MDET: block

MD
4
0.950060 1 0
-0.180171 2 1
-0.180171 2 1
-0.180171 2 1
DET 2 1 PR 2 1 3 1
DET 2 2 PR 2 1 3 1
DET 3 1 PR 2 1 4 1
DET 3 2 PR 2 1 4 1
DET 4 1 PR 2 1 5 1
DET 4 2 PR 2 1 5 1
END MDET

gwfn.data.gz

(1.32 KiB) Downloaded 42 times

1. Renormalize to c₁ = 1 for convenience. CASSCF then gives c₂ = −0.1896.
2. Pick a grid of |c₂| from 0 to 0.25 and freeze c₂ at each point.
3. At each point, optimize Jastrow + backflow by emin, as if the backflow were not quite right. This gives a one-parameter family of optimized nodes of controlled quality.
4. Compute the DMC energy at each point, extrapolated to dt → 0 at fixed N_w.
5. Evaluate Eq. (18) of Mitas & Annaberdiyev on the same wave functions.

Expected result: E^nda(c₂) should follow E_DMC(c₂) and have its minimum at the same c₂. The DMC optimum is near |c₂| ≈ 0.17, about 2.4 mHa below the HF node at c₂ = 0. If the two minima do not coincide, the displacement, measured in E_DMC, is the cost of using Eq. (18) instead of DMC to choose the node.

Variance of Eq. (18). It comes from two terms:
- Volume term ⟨V − V_Φ⟩. This is the local energy of Φ taken as a bosonic trial function. Its variance depends only on how close Φ is to the bosonic ground state, and it vanishes when Φ is exact.
- Surface integral. It is estimated from the configurations inside a tube of width ε around the node, sampled from Φ|Ψ|. Only a fraction of order ε² of the samples lands in the tube, so the relative error goes as √(1.5/n_tube) rather than 1/√N. A narrow tube lowers the bias but raises the noise. The measure Φ|Ψ| sets both.

A naive implementation is enough for an estimator, not for an optimizer. By a naive implementation I mean Φ = J_B·Π exp(−Z r_iI) and the tube estimate of the surface integral.
- On Be the tube holds ~3·10⁻⁴ of the sample, so the relative error of the surface term is ~56/√N, and 1% already needs ~3·10⁷ configurations.
- A backflow optimizer needs that accuracy on every step, for tens to hundreds of parameters, at the 0.1 mHa level of E_FN. On top of that, the estimator is only first-order accurate for a trial Ψ, so its minimum is displaced.

As an estimator, though, the same machinery also measures how hard the calculation itself is. For a given system it tells us:
- how much of the sample lands in the tube;
- how large the volume variance of Φ is;
- how the bias depends on ε and on Φ.

From these numbers follows the sample size needed to resolve a given difference in E_FN. Comparing that size with the cost of DMC tells us in advance whether the method is practical for that system at all.

[Mike Towler]

Vladimir, sounds like a good plan and the variance arithmetic is right; one thing. The ε² comes from the measure, Φ|Ψ| vanishing linearly at a simple node, not from the node itself, so it is a sampling choice rather than a property of the problem: a distribution that doesn'tt vanish there, |Ψ|^p with p below one, reweighted, makes the tube fraction of order ε and buys back a factor of order √(1/ε) in the surface term's error at the same N. Good luck with it; I'll look for the two minima.

[Vladimir_Konjkov]

Mike,

thanks for the |Psi|^p point, noted for the surface term. Before getting there I ran the DMC side of the protocol, and the first pass did not work out.

The system is Be with a 4-determinant CASSCF(2,4) wave function. The three excited determinants enter with a fixed coefficient c_i/c_0 = -C, and C is scanned over 0.00 ... 0.25. For each C the Jastrow and the backflow are optimized by emin. Every C is run at four time steps with dtdmc*nstep fixed (the same projection time on every rung), dmc_target_weight 1024, limdmc 4.

Attached are Be_4det.dat (columns C, dtdmc, E, error, measured stat inefficiency from POPSTATS) and energy.sh, which does the dt -> 0 extrapolation.

energy.png (50.92 KiB) Viewed 151 times

The best fit is linear in dt with one slope averaged over the whole scan. The per-C linear slopes agree within their errors, and the sqrt(dt) form buys nothing in chi2: E0 and c1 correlate at -0.98 to -0.99, which only inflates the error on E0. The global fit gives c2 = -0.00400 +/- 0.00042 and chi2/dof = 1.29 (25 dof):

Code: Select all

 C      E(dt->0)          C      E(dt->0)
0.00  -14.665080(31)    0.12  -14.667150(23)
0.01  -14.665205(31)    0.15  -14.667302(11)
0.02  -14.665565(39)    0.20  -14.667295(11)
0.05  -14.666247(22)    0.25  -14.666977(18)
0.10  -14.667010(26)

C = 0.15 and 0.20 form a plateau about 0.06 mHa above the exact -14.66736; they differ by 0.007 +/- 0.016 mHa.

The failure is in the points with bad efficiency. At C <= 0.12, eleven rungs have a measured stat inefficiency far from 1 (from -9.0 to +6.3), and their error bars are up to 0.94 mHa instead of the 0.03-0.05 mHa of their neighbours. For C = 0.15-0.25 every rung is clean. In my opinion this comes from a badly optimized backflow, not from DMC itself. Every C has its own emin backflow, about 130 parameters, and CASINO differentiates backflow parameters numerically, one full local energy per parameter. A backflow that the optimizer left off its minimum puts the node in a worse place, and that is exactly where the near-node excursions you described would come from.

So I am now writing analytical parameter derivatives of the backflow for CASINO. Only the Jastrow has them now; Slater and backflow parameters have has_aderiv = .false. What it takes:

1. Third derivatives of the orbitals: gaussians (gm4_bf in gauss_mol_bf) and STOs (the orbtderivs path of sto_orb_eval in stowfdet), a d3 helper in numerical.f90, and wfdet_has_tderivs in wfdet_basis to gate them.
2. Slater: the third derivatives of the determinant contracted with the backflow Jacobian (get_Tarray in slater.f90). The full tensor is never built: each term of d3(det)/det = L3 + 3 L2 L1 + L1 L1 L1 is contracted with D_ii as it is formed. This is the "tressian" algorithm described here: https://casinoqmc.readthedocs.io/en/lat ... l#tressian
3. wfn_aderiv_slater: the derivatives of ln Psi_S, of its gradient and of its Laplacian with respect to one backflow parameter, from the first, second and third derivatives with respect to the backflow coordinates and from the derivatives of X, dX and d2X with respect to the parameter. wfn_aderiv in wfn_utils dispatches to it.
4. pbackflow: has_aderiv for the expansion coefficients (the cutoffs stay numerical), with a FORCE_NUMDERIVS switch for A/B checks. The derivatives of X, dX and d2X come in one pass over the particles for every coefficient of eta, mu and Phi/Theta, and are cached per configuration. For a given parameter they are then contracted with the constrained change of the coefficients (unit_bf_param). The displacement and the cusp constraints are both linear in the coefficients, so this is exact. My first version ran the whole backflow_r2x again for every parameter, and on Be it gained only 1.5x.

Limits: real, non-pairing, 3D wave functions and polynomial backflow. Pi and Omega still go the per-parameter way; geminals and Pfaffians are not covered.

Once it is validated against FORCE_NUMDERIVS I will re-optimize the small-C backflows and redo the bad rungs, then come back to the two nodes.

Vladimir

data.tgz

(2.37 KiB) Not downloaded yet

[Mike Towler]

Vladimir, a plateau at 0.06 mHa above exact with a linear slope of −0.004 au⁻¹ is a nice result. One thing before you rebuild the backflow on the strength of the small-C rungs: POPSTATS's measured inefficiency isn't sign-definite. Its denominator is a difference of two nearly equal variance estimates taken under two weightings, one term of which can go either way, so it returns negative and wildly positive values with every walker weight positive, and the guard fires only at +2 or above. I have just been through that maaths in the code for the paper. Your bad rungs' error bars are real; the inefficiency number next to them doesn't say anything about the node though. Analytic backflow parameter derivatives would be very welcome..