{'fraction_non_alphanumeric': 0.035438884331419195, 'fraction_numerical': 0.017596390484003282, 'mean_word_length': 4.2920555676145, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective noise reduction techniques for disconnected loops in Lattice QCD

Many Lattice QCD observables of phenomenological interest include so-called all-to-all propagators. The computation of these requires prohibitively large computational resources, unless they are estimated stochastically. This is usually done. However, the computational demand can often be further reduced by one order of magnitude by implementing sophisticated unbiased noise reduction techniques. We combine both well known and novel methods that can be applied to a wide range of problems. We concentrate on calculating disconnected contributions to nucleon structure functions, as one realistic benchmark example. In particular we determine the strangeness contributions to the nucleon,, and to the spin of the nucleon, Delta s.


Introduction
Advanced Lattice QCD calculations often require the evaluation of diagrams with disconnected quark lines. Important examples are properties of flavour singlet mesons [1], QCD spectroscopy including multiquark and scattering states [2], the determination of hadronic scattering lengths [3], and of isosinglet contributions to hadronic form factors and structure functions [4,5]. Moreover, statistical errors of some standard observables like meson masses and electroweak decay constants can be reduced by averaging the sources over the lattice volume. In all these cases the standard point-to-all propagators that are obtained by calculating twelve (colour times spin) columns of the inverse lattice Dirac operator are not sufficient. Instead, all-to-all or timeslice-toall propagators need to be computed.
Lattices typically contain a number of sites ranging from V ≈ 10 6 up to V ≈ 10 9 . Inverting a lattice Dirac operator by conventional means would increase the effort in terms of computer memory and operations spent by this factor, relative to the cost of obtaining a single point-to-all propagator. An alternative approach to the problem consists of calculating an unbiased stochastic estimate of the propagator, replacing a factor 12V in effort by a number N of random sources, where in certain situations N can be as small as 10. Such estimation is permissible since, provided it is unbiased, it will only add to the statistical errors. These are present in any case, due to the path integral evaluation of expectation values of observables by means of a Markov Monte Carlo simulation (for an introduction into Lattice QCD simulations see e.g. [6]).
In general, the additional error introduced by this procedure will, for sufficiently large numbers of estimates N, decrease in proportion to 1/ √ N. Ideally, N should be chosen such that this additional error is smaller than the statistical error induced by the Monte Carlo time series. The optimal number will depend on the observable in question and on the prescription employed to calculate the stochastic estimate. In this article we discuss different improvement methods aimed at decreasing the prefactor of the asymptotic 1/ √ N behaviour. In particular, we will introduce a new such method which we call the truncated solver method (TSM). We then test combinations of different improved stochastic estimator methods in a realistic Lattice QCD simulation. As our benchmark examples we choose to calculate the strangeness contribution to the spin of the nucleon ∆s as well as the scalar strangeness content N|ss|N .
The spin of the nucleon can be factorized into a quark spin contribution ∆Σ, a quark angular momentum contribution L q and a gluonic contribution (spin and angular momentum) ∆G: In the naïve nonrelativistic SU(6) quark model, ∆Σ = 1, with vanishing angular momentum and gluon contributions. In this case sea quark contributions will be absent too and therefore there will be no strangeness contribution ∆s in the factorization, where in our notation ∆q contains both, the spin of the quarks q and of the antiquarksq. Experimentally, ∆s is usually obtained by integrating the strangeness contribution to the spin structure function g 1 over the momentum fraction x. The integral over the range in which data exists (x 0.004) usually agrees with zero. For instance a recent Hermes measurement in the region x ≥ 0.02 yields [7] ∆s = 0.037 (19) (27). This means that non-zero results rely on extrapolations into the experimentally unprobed region of very small x and are model dependent [8,9]. The standard Hermes analysis [10] yields ∆s = −0.085(13)(8)(9) at a renormalization scale µ 2 = 5 GeV 2 in the MS scheme. Our results below suggest, ∆s = −0.017 (21), unless there are large effects in the chiral extrapolation from a pseudoscalar mass m PS ≈ 450 MeV to the physical point.
The scalar strangeness density is not directly accessible in experiment but plays a rôle in models of nuclear structure. It is also of phenomenological interest since, assuming that heavy flavours are strongly suppressed, the dominant coupling of the Higgs particle to the nucleon will be accompanied by this scalar matrix element.
The outline of this article is as follows: in Sec. 2 we introduce our notation and explain the stochastic estimator and noise reduction methods applied. In Sec. 3 we detail our lattice setup and employ these techniques to calculate disconnected quark loop contributions. Finally, in Sec. 4 we present our results on the axial and scalar nucleon strangeness matrix elements, before we conclude. The systematic tuning of the parameters used in the truncated solver method is described in detail in Appendix A. The TSM is quite generally applicable. This is also demonstrated in the appendix where the conjugate gradient (CG) solver is replaced by the popular stabilized biconjugate gradient (BiCGstab2) solver [11]. Preliminary results appeared in [12] and [13].

Noise reduction techniques
As outlined above we require propagators M −1 connecting arbitrary pairs of lattice points, where M denotes a lattice Dirac operator. In our specific case we employ the Wilson quark action [14,6]. Since the propagators have 12V×12V components, where V denotes the number of lattice points, direct evaluation would be prohibitively expensive, both in terms of computer time and of memory. However, we encounter statistical errors anyway from the importance sampling of path integral expectation values. Hence it is sufficient to aim at an (unbiased) estimate, which can be obtained by stochastic methods [15]. For this purpose we introduce the following notation, N denotes the number of "stochastic estimates". Let |η i , i = 1, . . . , N be random vectors with the properties, These requirements are for instance met by complex 2 noise where the 12V components are numbers e iφ , with the uncorrelated random phases φ ∈ {±π/4, ±3π/4}. In [16,17,18] it has been demonstrated that real and complex 2 noise reduces the variance, relative to other choices such as Gaussian or double hump noise. In our experience, similarly small variances can also be obtained with N , U(1) or even with SU(3) noisy sources [19], indicating that the only important property is an equal modulus of the components. In the present study we employ complex 2 noise, where the components of the random vectors |η i run over the spacetime volume, spin and colour. We define the Hermitian lattice Dirac operator Q = γ 5 M. If we solve the linear systems, for |s i , then we can substitute, see Eq. (5), The difference between the approximation of Eq. (7) above and the exact result is unbiased and reduces like 1/ √ N. The sparse linear system of Eq. (6) can for example be solved by means of the conjugate gradient (CG) or the stabilized biconjugate gradient (BiCGstab2) algorithms. Note that in the CG case we can actually use even/odd preconditioning by employing the operator M † M = Q 2 . We then obtain |s i by multiplying the result with Q.
The stochastic estimate approach reduces the problem from O(12V × 12V) to O(N × 12V), in terms of memory and computer time. The stochastic error will remain roughly constant if N is scaled like √ Va 2 with the lattice volume, where a denotes the lattice spacing. In order to limit the computational effort, N should not be chosen overly large. However, in the end the stochastic noise should not be the dominant source of statistical error. In general, the optimal balance between the stochastic sampling on single configurations and the Monte Carlo importance sampling of gauge configurations will depend on the observable in question and on the methods used. 3 Due to the difference between |s η| and Q −1 above, any fermionic observable A can only be estimated up to a stochastic error ∆ A,stoch = O(1/ √ N) on a given configuration. We define the configuration average · c over n uncorrelated configurations and normalize this appropriately: For large N and n this will scale like σ 2 A,stoch ∝ (Nn) −1 . We also define the total error σ A,tot as the variation of the obtained estimates of A over gauge configurations. At fixed N this will be proportional to 1/ √ n. Obviously, this total error is limited by the stochastic error 1 : If σ 2 A,stoch ≃ σ 2 A,tot then it is worthwhile to improve the quality of the estimates while if σ 2 A,stoch ≪ σ 2 A,tot then precision can only be gained by increasing the number of configurations n, possibly reducing N to save computer time since the n −1 scaling cancels from the above inequality.
These methods are explained below.

Partitioning
We decompose R = V ⊗ colour ⊗ spin into m subspaces R j : R = ⊕ m j=1 R j . One can then set the source vectors |η i | j to zero, outside of the subspace R j . We label the corresponding solutions as |s i | j . The solution for the all-to-all propagator is then given by the sum, This procedure results in an m-fold increase in the total number of solver applications. The term ½ − |η η| multiplying Q −1 in Eq. (7) only has off-diagonal entries, all of similar sizes. Therefore the off-diagonal entries of Q −1 will determine the stochastic error of the particular observables in question. Q −1 will exponentially decay with the spacetime distance between source and sink: spacetime components in the neighbourhood of the sink position will yield the leading contributions to the stochastic variance. Likewise the noise for spin components that are strongly coupled to each other by Q −1 , multiplied by the relevant Γ and derivative structures, will dominate. Ideally, partitioning will black out the largest contributions to the stochastic noise. If the achieved reduction exceeds the 1/ √ m factor, then the computational overhead is justified. This overhead 1 For an exact calculation of A (σ A,stoch = 0), a Monte Carlo error σ A,stat = σ A,tot ∝ 1/ √ n can be introduced. In the limit of large N and n, the factorization σ 2 A,stat = σ 2 A,tot − σ 2 A,stoch holds, see e.g. [20]. For our considerations there is no need to introduce this quantity.
can be somewhat reduced at the expense of memory by sophisticated preconditioning techniques and/or aggressive deflation [27,28,29,30,31,32]: the constant set-up costs of such techniques become more affordable with a larger number of right hand sides. We experimented with various patterns and combinations of colour, spin and time partitioning, see e.g. [33]. Partitioning as a stand-alone solution works very well in many situations [21,34,20,22]. However, such gains are mostly eliminated once partitioning is combined with the other tricks. Time partitioning is a notable exception: in situations where only timeslice-to-all propagators are required there is no increase in the number of solves but the variance is still reduced. In the calculation of disconnected contributions to the nucleon structure it also turns out to be useful to generate the current insertion at more than one timeslice, e.g. to average the correlation with a nucleon propagating from a given source in the forward direction and the backward propagating one (which comes for free). In this case one can seed the random sources at two (or more) well separated timeslices and still gain from partitioning, without any associated overhead.

The hopping parameter expansion technique
The stochastic noise from terms that are close to the diagonal is accompanied by larger amplitudes than terms that are far off the diagonal, see Eq. (7) and the discussion in Sec. 2.1 above. Hence the cancellation of near-diagonal noise requires a comparatively larger number of estimates. The HPE aims at eliminating some of the near-diagonal noise contributions by exploiting the ultra-locality of the action. Thus it cannot be generalized for instance to the Neuberger action [35].
We rewrite the fermionic matrix as, The HPE is based on the observation that one can expand, where k ≥ 1. For distances between source and sink consisting of more than k links, the first term on the right hand side does not contribute since D only connects nearest spacetime neighbours. Therefore, M −1 xy = [(κD) k M −1 ] xy for sufficiently large source and sink separations. However, at finite N, where the variance of the latter estimate of M −1 xy is reduced. This was for instance exploited in [25]. In some cases additional powers of D can be gained due to the Γ structure of the creation and annihilation operators.
Here we study closed loops, i.e. x = y. Obviously, only even powers of D contribute to where the first term can be computed trivially. In the other cases, For the Wilson action, k = 8 for Γ ∈ {γ µ γ 5 , γ 5 } and k = 4 otherwise. For the Sheikholeslami-Wohlert action [36], k = 2. The lowest non-vanishing terms have been calculated analytically [23] and can be computed and corrected for exactly (unbiased noise subtraction) [23,37,5].
We limit ourselves to the highest vanishing order (k = 4 or k = 8 for the Wilson action), where such subtraction is not necessary. In this case the computational overhead of the HPE is small, such that this substitution can easily be combined with the TSM that we explain below.

The truncated solver method
It has been noticed long ago [38] that solvers typically converge to the correct result within an accuracy of the size of the stochastic error after a relatively small number of iterations. Practitioners have therefore sometimes relaxed the requirement on the residual when solving for stochastic sources. This is not unproblematic since it introduces a systematic bias that will be invisible on one configuration but might very well affect the result obtained on a sample of, say, 200 configurations. We label the result obtained after n t solver iterations by |s i (n t ) , where omitting the subscript means convergence within numerical accuracy. We can now factorize: The above equation is an exact linear decomposition of Q −1 and the algorithm used to calculate both parts is well defined. Due to these properties and the fact that the |η i with i ≤ N 1 are uncorrelated with those for i > N 1 , the resulting estimate is unbiased. Ideally, one will generate a large number N 1 ≫ N 2 of relatively cheap estimates at small n t and then remove the bias by correcting with a small number N 2 of expensive solutions to machine precision. This method can easily be combined with all the other methods described here. It is possible to tune the two parameters N 1 /N 2 and n t that enter the algorithm to the point of minimal variance at fixed computer time in a relatively inexpensive and straightforward way. This is discussed in Appendix A.
In some sense the underlying philosophy of TSM is similar to estimating the cheap first term within the HPE Eq. (12) with many random sources and the expensive second term with only a few sources. However, iterative solvers like the CG converge faster than the HPE. Moreover, TSM is applicable to any fermion action, not only to ultra-local ones, and efficient for any quark mass. TSM can be combined with HPE, see Sec. 2.5 below.
In Appendix A we also demonstrate that the TSM can be generalized to other solvers, by employing BiCGstab2. However, the smooth convergence of CG, that is relatively independent of the random source and gauge configuration, is a clear advantage when it comes to deciding on a TSM parameter set. Moreover, in the CG case the combination with the truncated eigenmode acceleration bears less computational overhead and is compatible with even-odd preconditioning (see below). We emphasize however that the converged solutions |s i within Eq. (13) can be obtained by using any efficient solver. In particular, there is no need to employ the same solver as that used for obtaining the truncated solutions |s i (n t ) . The truncated solver of course needs to be the same for i > N 1 as that used for i ≤ N 1 .

The truncated eigenmode acceleration
We define the n e smallest (real) eigenvalues q i of Q and the corresponding orthonormal eigenvectors |u i , i = 1, . . . , n: These can for instance be calculated by means of the parallel implicitly restarted Arnoldi method (IRAM) with Chebychev acceleration [39] or by Ritz methods [40]. 6 We can now approximate [39], However, this estimate is biased. We define the projection operator È n e , onto the complement of the space spanned by these n e eigenvectors, Projecting the stochastic sources onto the orthogonal complement, yields the new unbiased estimate [25], where a left projection of |s i ⊥ is not necessary since [Q, È n e ] = 0 and thus È n e Q È n e = Q È n e . The above estimate will usually have a reduced variance. Note that in the literature exact point-to-all propagators have also been combined with a truncated low eigenmode all-to-all calculation to achieve smaller gauge errors (low mode averaging) [41,42].
A nice side effect of TEA lies in the acceleration of the solver, by deflation [43,27,29,31,44,32]. The condition number of the projected operator Q È n e and therefore the number of solver iterations are reduced. Within some algorithms like the CG the fact that the projector commutes with the operator Q guarantees that the Krylov subspace remains confined within the orthogonal complement. So in these cases no further intermediate projections are necessary to fully exploit the potential of deflation. Note that while the number of undeflated solver iterations increases like 1/m 2 PS at small quark masses, the efficiency of the eigenvalue calculation remains the same. Unfortunately, at small m PS , one would like to increase the linear lattice extent L ∝ 1/m PS . In this case, the rank of the deflation space in the worst case will increase like n e ∝ L 4 ∝ Va 4 . Depending on the volume required this may or may not become a serious problem. The volume scaling of stochastic methods is somewhat more favourable: the number of estimates needs to be adjusted at most like √ Va 2 (this can be somewhat reduced by partitioning). From these considerations it becomes evident that the optimal choices of n e and N strongly depend on the volume, the quark mass and the lattice spacing. The same holds for the algorithm where a deflated CG can outperform a (deflated) BiCGstab2, in particular when combined with the truncated solver method.
Obviously, it is also possible to decompose M, rather than the Hermitian operator Q = γ 5 M of Eq. (15), into eigenmodes [45]: M|r i = λ i |r i . However, in this case we end up with a biorthonormal system ℓ i |r j = δ i j , where the left eigenvectors |ℓ i will differ from the right eigenvectors: ℓ i |M = ℓ i |λ i . One can now decompose, It can easily be seen that due to the property M † = γ 5 Mγ 5 ,λ i := λ * i is an eigenvalue whenever λ i is an eigenvalue, with a left eigenvector l i | = r i |γ 5 and a right eigenvector |r i = γ 5 |ℓ i . The advantage of this decomposition is that the eigenvectors are independent of the quark mass and the eigenvalues at different κ values are trivially related to each other. This follows from the structure Eq. (11). In this article we will not pursue this alternative eigenmode decomposition any further. However, the decomposition of M −1 might converge better in some channels than that of Q −1 and vice versa. While this biorthogonal approach is incompatible with deflating the CG solver, it would be the natural starting point for BiCG solvers.

Combining the methods
Partitioning can trivially be combined with the other three methods. However some notes are in place with respect to combinations of HPE, TSM and TEA. Within the HPE a leftmultiplication of the estimate E[Q −1 ] with D k is essential since we have defined M −1 = Q −1 γ 5 . A right-multiplication would involve γ 5 Dγ 5 = D † . The easiest way of implementing HPE is to multiply the solution vectors with (κD) k , prior to any application of smearing functions or contractions but after TEA and TSM.
Within Eq. (18) only the projected source vectors appear. Hence after the projection, Eq. (17), the original noise vectors |η i can be discarded. The overhead from the projection can be significant when combined with the TSM, Eq. (13), where the large number of estimates N 1 implies a large number of projections and the small number of solver iterations n t indicates that not much computer time will be spent between these projections. The projection overhead can be reduced by restricting the computation of the inner products to the partitioning subspace: many components of |η i might have been set to zero if the partitioning method was used. This saving cannot be obtained at the sink. Fortunately, [È n ev , Q] = 0, such that the solutions |s i ⊥ remain within the orthogonal complement and no second projection is necessary, provided the residual of the solver is chosen sufficiently small! Obviously the residual will not be small for the truncated solutions, which therefore one would wish to project back into the orthogonal complement (before application of the D operator within the HPE). Fortunately, the (even-odd preconditioned) CG algorithm never leaves the orthogonal subspace: in this case, such an additional projection is never necessary, even if the solver is not run to convergence. Neither is such a left projection strictly required for solvers without this property, as long as the precision of the part that is run to convergence is sufficiently large. In this case, the N 1 biased estimates will pick up some unwanted low mode contributions that will however be corrected for by the N 2 estimates of the bias. This might or might not increase the stochastic errors.
For completeness we mention that domain decomposition techniques have been suggested in the literature [46,47]. These can easily be combined with TSM and TEA too, as can the so-called one-end-trick [48,49,50]. We remark that if a bosonic representation is employed [46,51], estimating M † M = Q 2 instead of Q, then the TSM can in principle be substituted by a quasi heatbath strategy [52].

Lattice setup
We study combinations of the variance reduction techniques outlined above, using configurations provided by the Wuppertal group: these are n f = 2 + 1 dynamical configurations with V = 16 3 × 32 lattice points, generated using a Symanzik improved gauge action and a stout-link improved (rooted) staggered fermion action. The lattice spacing is fairly coarse, a −1 ≈ 1.55 GeV, while the spatial extent is around 2 fm. Further details can be found in [53].

Results for the disconnected loop
We wish to calculate nucleon structure observables. For this purpose a nucleon will be created at some initial time t 0 and destroyed at a later time t f ≫ t 0 , with a current inserted at some intermediate time t. This is illustrated in Fig. 1. The disconnected loop Tr (M −1 Γ) within the right diagram of the figure will be calculated with stochastic all-to-all techniques and can subsequently be combined with the nucleon two-point function, calculated in the standard point-to-all way, on a configuration by configuration basis.
A preliminary picture of how well the variance reduction techniques work can be obtained by studying the zero-momentum projected disconnected loop We investigate the reduction in computer time of optimized stochastic estimates, relative to the situation without the use of any improvement techniques (except for time partitioning: all disconnected loops are calculated on one timeslice only). We state all costs in terms of the average real computer time required on a Pentium 4 PC for one solver application (unimproved estimate), where we account for all overheads of the improvement methods. We also ran the parallelized code on IBM p6 and Opteron clusters as well as on BlueGene/L and BlueGene/P systems. The overall efficiencies of matrix-vector multiplies and global sums depend on the architecture, however, the gain factors defined below are not significantly affected.
We define the gain as the ratio, taken at fixed real computer time. var denotes the variance between stochastic estimates on a given gauge configuration, E 0 and E imp stand for the unimproved and improved estimated, respectively. We start by investigating the TSM at the heaviest quark mass, κ loop = 0.166. Within TSM the combination Tr (M −1 Γ) is obtained as an average over N 1 truncated solutions |s i (n t ) : η i |γ 5 Γ|s i (n t ) . This value is then corrected by N 2 ≪ N 1 estimates of the resulting bias, see Eq. (13). The faster the truncated estimate as a function of n t approaches the estimate obtained at full convergence, the better the method will work. In this study we only consider currents of local quark bilinears, where we use the Γ conventions of [54]. We observe very satisfactory convergence rates for all the 16 possible Γ structures. In Fig. 2 we illustrate this for the even/odd preconditioned CG solver for the worst case (Γ = ½) and for the best case (Γ = γ 3 γ 5 ). Note, however, that the unimproved estimate of Γ = ½ is already very precise to start with. Convergence at this κ value is reached after n conv ≈ 480 CG iterations.
Applying TSM involves making choices for n t and for N 1 /N 2 . We detail our optimization procedure in Appendix A. This systematic tuning results in similar amounts of computer time being spent on the truncated estimate as on estimating the bias. In particular we find N 1 /N 2 ≈ n conv /n t . Performing this optimization on a single configuration appears sufficient since the variance is not 10  much affected if n t and N 1 /N 2 are changed by 20 %. Fluctuations between configurations turn out to be smaller than this value. In Table 1 we display the optimized parameters at κ loop = 0.166 for a representative choice of Γ structures (scalar, vector, tensor, pseudoscalar and axial vector), together with the approximate fixed cost gains, Eq. (21). These factors vary between 5 and 10. In a real production run one would not wish to generate new sets of optimized estimates for each Γ structure or observable one is interested in. Fortunately, with the still tolerable exception of Γ = ½, n t and N 1 /N 2 exhibit only a mild Γ dependence such that similar overall gains can still be achieved with just one parameter set.
TSM can trivially be combined with the hopping parameter expansion. At κ loop = 0.166 the gains from applying a stand-alone HPE are almost as big as those from the TSM. However, the HPE does not work well for Γ = ½, see Fig. 3. Both methods aim at removing noisy short range contributions from the estimates. Clearly, the gain from combining the two methods will be smaller than the product of the two isolated gains. However, the separation of long and short range is organized differently. The HPE only works for ultra-local actions and explicitly removes terms up to a particular lattice hopping radius. It will be less convergent and hence less effective at very light quark masses. The TSM is more generally applicable and separates short range modes from long range ones, but this is done within the Krylov space of the solver. When combining HPE with TSM, the reduced variances of the truncated estimates result in larger optimal n t values. The constant computational overhead per solve of the additional D applications also pushes n t gain TSM TSM+HPE m PS ½ γ 3 γ 1 γ 2 γ 5 γ 3 γ 5 ½ γ 3 γ 1 γ 2 γ 5 γ 3 γ 5 600 MeV 5 5    towards larger values. Table 1 demonstrates that the combined gain of these two methods at κ loop = 0.166 can be as large as a factor of 30.
Having tested the method at κ loop = 0.166, where m PS ≈ 600 MeV, we also use it to calculate the disconnected loop at κ loop = 0.1675 and at κ loop = 0.1684, corresponding to pseudoscalars masses of approximately 450 MeV and 300 MeV, respectively. The resulting gains are displayed in Table 2. While the TSM performance is very independent of the quark mass, the HPE becomes less effective at lighter masses, in agreement with our expectation. Still, even at the lightest quark mass, the combined gains range from factors of 6 (for Γ = ½) up to 19 (for Γ = γ 3 γ 5 ).
At light quark masses where the HPE becomes less effective, the low eigenmode contributions to hadronic observables might become more dominant [39,41,42,25]. To study this effect, we deflate the lowest eigenmodes at the lightest pseudoscalar mass. As expected, this accelerates the solver [43,26,27,28,29,30,31,44,32]. We display a typical residual, as a function of the number of CG iterations, without deflation and deflating the lowest 5, 10, 20 and 50 modes of the Hermitian Wilson-Dirac operator in Fig. 4.
The optimized TSM parameters at κ loop = 0.1684 for combining TSM with HPE as well as 12   for combining TSM with HPE and the TEA of the lowest 20 eigenmodes are shown in Table 3. The faster rate of convergence of the deflated solver leads to smaller n t values and therefore to larger N 1 /N 2 ratios when TSM is combined with TEA. When normalized to the real cost of generating 300 unimproved stochastic estimates, all gain factors increase by more than a factor of two, and this in spite of one quarter of the time being spent on generating and projecting onto the eigenvectors. However, this factor can fully be attributed to the acceleration of the solver: while there exist observables where stochastic errors decrease when applying the TEA [39,41,25], for the quark loops that we investigate here this is not the case, at least not at pseudoscalar masses above 300 MeV. We observe a break-even between TSM+HPE+TEA(È 20 ) and TSM+HPE when matching the cost of approximately 100 estimates, as can be seen in the last row of the table. In lattice simulations we will encounter gauge errors from the Monte Carlo time series, in addition to the errors from the stochastic estimates on single configurations, discussed here. This interplay is studied below.

Configuration averages
When calculating averages over configurations, the stochastic errors σ stoch , that are defined in Eq. (8), should be smaller than say half of the final Monte Carlo errors σ tot . However, making σ stoch unnecessarily small can be a waste of computer time that would be more wisely spent on analysing more configurations or realizing additional source positions. We start to study this balance of errors by computing expectation values over configurations for the disconnected loops. In Sec. 4 below we will then calculate the three point functions that are of physical interest.
We remark that the (zero momentum projected) gauge average Re Tr (M −1 ½) c will have a non-trivial value, while Re Tr (M −1 γ µ γ 5 ) c = 0. We first assess the impact of low mode deflation. In Fig. 5  and deflating 20 modes at κ = 0.1684 results in a two third reduction (red versus blue symbols of the right figure). We observe variations of n conv between and within configurations that increase with decreasing quark mass, as has been investigated quantitatively in [57]. Deflating the lowest modes vastly reduces the variance in n conv . Clearly, the TSM parameter tuning benefits from this stabilization.
We now average our improved estimates over n configurations, where we employ TSM, HPE and, at the lightest κ value, TEA with a projection onto the 20 lowest modes. In addition, we calculate up to 100 unimproved estimates on each configuration. The results as functions of the real computer time spent, in units of the cost of one unimproved estimate, i.e. of one solve, are displayed in Table 4, for Γ = ½ and Γ = γ j γ 5 , where we average over the three possible j-directions. The total errors σ tot are displayed in brackets, followed by their respective lower bounds, as given by the stochastic errors σ stoch .
For Γ = ½, even at the cost equivalent of six standard estimates, the gauge errors still dominate over the improved stochastic errors. Therefore, there is little point in exceeding this number. On the given ensemble, the same error can be obtained using 25 to 50 standard estimates, suggesting a five-fold saving in computer time. However, this is somewhat misleading since, by just increasing the number of gauge configurations by a factor of 1.6, the same error can be obtained, employing six unimproved estimates. It should be noted that the error balance could in principle look differently, once the loop is correlated with a nucleon propagator, within a three point function. Moreover, the stochastic error will become more relevant on large volumes. At the lightest κ value, due to the overhead of setting up TEA, the cost of the improved estimates will always exceed the number 6. In fact, as can also be seen from Table 4, TEA only turns out to be a worthwhile enterprise for a cost equivalent bigger than 90. Hence, unless the eigenvectors have been generated anyway, for instance to enable low mode averaging [41] for the two point functions or for the calculation of other current insertions, TEA is best omitted for Γ = ½.
The picture for the axial current containing Γ = γ j γ 5 is different: here at the cost equivalent of 100 standard estimates, the stochastic error still accounts for one quarter of the total error and the gain of applying the improved method is in all cases more than four-fold. This advantage should increase further at larger volumes. Note that the reductions of the stochastic errors alone, for the scalar and axial cases, are consistent with the results of Table 2 that were obtained on a single configuration.
Based on these results, we decide to invest the cost equivalent of 100 even/odd preconditioned CG solves into the TSM/HPE/(TEA) estimation of the axial loop and of 6 CG solves into the estimation of the scalar loop. Note that for the calculation of a standard point-to-all baryonic two-point correlation function usually 12 such solves are required, and even more sources if a variational basis is used to optimize the creation operator. We display the resulting stochastic 14 errors ∆ stoch on single configurations at κ loop = 0.166 in Fig. 6. This graphical representation again makes it evident that even for this low-cost setting the stochastic errors are much smaller than the gauge noise. The larger error bars on the right of the figures correspond to unimproved ∆ stoch values, obtained at the same cost. After application of the improvement methods, the stochastic errors are so small that nothing extra can be gained from increasing the number of stochastic sources (and solves) beyond these moderate values.

Definition of the matrix elements
We now apply our methods to the calculation of observables of phenomenological interest, namely of ∆q and N|qq|N . The contribution to the nucleon spin ∆q is defined through the matrix element (in Minkowski space notation), where m N denotes the mass of the nucleon N and s µ its spin (s 2 µ = −1). ∆q and N|qq|N are extracted from the ratios of three-point functions to two-point functions (at zero momentum): For the scalar matrix element we use, Γ 2pt = Γ unpol := (1 + γ 4 )/2 and Γ loop = ½. For ∆q we calculate the difference between two polarizations: Γ 2pt = γ j γ 5 Γ unpol and Γ loop = γ j γ 5 , where we average over all three possible j-orientations. The spin projection operators along the jaxis read, P ↑↓ = 1 2 (½ ± iγ j γ 5 ), so that in this case, Γ 2pt = −i(P ↑ − P ↓ )Γ unpol , where we have traded a factor i against taking the real part, rather than the imaginary part, of the nominator in Eq. (23). The variance of the above expression is reduced by explicitly using the fact that Im Tr (M −1 ½) = Im Tr (M −1 γ j γ 5 ) = 0. The two-point function of the zero momentum projected proton with sink and source spinor indices α and β is given by, where we have set t 0 = 0. This can be constructed from standard point-to-all quark propagators [6]. Note that for q = u, d there are additional connected contributions R con , which we have not calculated. We combine the three κ loop values with κ 2pt = 0.166 and 0.1675. In the limit of large times, t f ≫ t ≫ 0, in the axial case, The prefactor two comes from taking the difference between two opposite polarizations rather than fixing one polarization. We note that this result will be in a lattice scheme and still needs to be multiplied by a renormalization factor of O(1), for a translation into the MS scheme. In the scalar case, N|qq|N is defined as the connected contribution only and can thus be obtained by subtracting 0|qq|0 = − x Re Tr (M −1 (x, t; x, t) c from Eq. (23). We will label the disconnected contributions to these two matrix elements as, ∆q dis and N|qq|N dis , respectively. In the case of the strange quark, ∆s = ∆s dis and N|ss|N = N|ss|N dis .

Results of the calculation
We calculate the disconnected loop on only one timeslice t = 3 ≈ 0.38 fm/a. Having the operator inserted close to the source (t 0 = 0) reduces the statistical errors but care must be taken 17 that contributions from excited states are suppressed so that Eq. (25) still holds. Following [25], we use Wuppertal smearing on top of APE smeared links at the source and sink for both the three point and two point functions. In Fig. 7 we display the effective mass, as a function of the time. This shows that with our choice of smearing parameters the excited state contributions are not significant at t = 3. Hence, in the symmetric situation, t f = 6, such contributions should be negligible. This is even more so since the statistical error will be much bigger than that of the two point function at t = 3. The consistency of this assumption can be checked by varying t f ≥ 4. In Fig. 8 we display finite time estimates of the matrix elements N|qq|N dis and ∆q dis in the lattice scheme, for different final times t f at κ loop = κ 2pt = 0.166. This value approximately corresponds to the mass of the strange quark. The stochastic estimates of the loops were generated using time partitioning (t = 3), TSM and HPE, at the cost of 6 even/odd preconditioned CG solves for the scalar and at the cost of 100 such solves for the axial current. Our analysis of the two point function above suggests that finite t systematics should be small, relative to the statistical errors at t f = 6. Indeed, all t f ≥ 4 data are consistent with a constant. We conservatively quote the t f = 6 values as our final results. We realized all six κ loop ∈ {0.166, 0.1675, 0.1684} and κ 2pt ∈ {0.166, 0.1675} combinations where the lightest κ loop value corresponds to m PS ≈ 300 MeV. The five combinations that are not shown result in the same general picture, with larger statistical errors. For the two point function we do not realize the lightest κ value since, without low mode averaging, this turns out to be very noisy.
In Table 5 we display improved and conventional stochastic estimates of the scalar matrix elements (obtained at t 0 = 0, t = 3 and t f = 6 ≈ 0.76 fm/a), for all our κ combinations. The fixed cost reductions, due to TSM and HPE, in the relative errors are small and do by far not match the gain factors that we obtained in Sec. 3 for the disconnected loops alone. The noise is dominated by taking the correlation with the two point function. For the precision of the disconnected loop to be matched by that of the two point function, the latter needs to be evaluated for multiple source points, eventually in combination with low mode averaging [41]. We note that one does not encounter any computational overhead in calculating the disconnected loop at more than one timeslice. As long as these are sufficiently separated, the stochastic errors will not increase 18   (20) 0.011 (38)  significantly.
In Table 6 we display the same information as in Table 5, for ∆q dis in the lattice scheme. In this case, at the cost equivalent of 100 estimates (90 estimates at κ loop = 0.1684), the gains of applying TSM and HPE in terms of computer time, relative to time partitioning alone, are about ten-fold. The reductions in error are close to those displayed in Table 4 for the disconnected loop alone. This channel is not yet limited by the accuracy of the two point function but of course also in this case statistics could be increased by averaging over multiple baryon sources and over forward as well as backward propagation.
In Fig. 9 we display our final t f = 6 results for the two matrix elements. In neither case do we observe any significant dependence on the valence quark mass, varying this from m PS ≈ 600 MeV down to 450 MeV, or on the loop quark mass, reducing m PS ≈ 600 MeV (strange quark mass) to m PS ≈ 300 MeV. We find ∆s = −0.020(11) at the heavier proton mass and ∆s = −0.017(21) at the lighter mass value, while the scalar matrix elements are somewhat larger than one. Note that the lattice results presented here are unrenormalized. However, based on perturbative two loop results [58], albeit with different sea quark and gluon actions, we would not expect ∆s to change by more than a factor of 0.7 after a translation into the MS scheme.

Conclusions and outlook
A growing number of Lattice QCD applications requires all-to-all propagators. These are most efficiently estimated stochastically. We presented the novel truncated solver method (TSM) that typically reduces the computer time to achieve a given stochastic variance by an order of magnitude, without introducing a bias. The gain factors of this method for different observables were demonstrated not to vary when changing the quark mass by a factor of four (600 MeV ≥ m PS ≥ 300 MeV). The TSM is easy to implement and can be used for any fermionic action. The combination of TSM with different variance reduction techniques is straight forward and this reduces the stochastic variance even further. We studied in detail combinations of TSM, partitioning [18], the hopping parameter expansion [23] and low eigenmode deflation [26,25].
In realistic Lattice QCD simulations there are intrinsic errors from the Monte Carlo time series, in addition to the errors introduced by the stochastic estimation of the inverse of the fermionic matrix on individual configurations. We studied the interplay between these gauge and stochastic noises. After reducing the stochastic variance by a combination of methods, in our calculation of disconnected contributions to the nucleon structure, the gauge errors became the dominant sources of uncertainty. This means that we can afford larger stochastic errors and for instance increase the lattice volume, without having to increase the number of random sources. For instance, on our 2 fm lattices, even at a cost of only 6 CG solves, the stochastic error of the scalar matrix element N|qq|N dis is completely over-shadowed by its gauge error: in certain situations one can overdo the reduction of the stochastic noise. In this particular case the total error can more efficiently be reduced by increasing the number of nucleon sources on each configuration than by increasing the number of (improved) estimates. Also the determination of ∆s will benefit from this. At present we are pursuing such an approach.
Our result on the strangeness contribution to the nucleon spin ∆s ≈ 0 is in agreement with the other recent direct calculation of this quantity [59]. However, we disagree with earlier, less precise studies that employed a summation method over t [60,61,62]. These suggested a value ∆s ≈ −0.12. In the case of the scalar matrix element our errors are still too large to state a meaningful number, in particular since chiral and infinite volume behaviours need to be studied. The techniques developed here are used by us in an ongoing study [63] at smaller lattice spacing, quark masses and larger volumes with Sheikholeslami-Wohlert sea and valence quarks.

Appendix A. TSM parameter tuning for the CG and BiCGstab2 solvers
The truncated solver method of Sec. 2.3 depends on two parameters, namely on the number of iterations for the inexact solves, n t , and on the ratio N 1 /N 2 of the number of inexact estimates over the number of estimates of the correction term, see Eq. (13). These parameters need to be fixed, ideally, so as to minimize the variance of the estimate of the disconnected loop, E[Tr (M −1 Γ)], at fixed cost. The estimates are uncorrelated and for N 1 , N 2 ≫ 1 the variance factorizes into, var[Tr (M −1 Γ)] = var η|γ 5 Γ|s (n t ) N 1 + var η|γ 5 Γ( |s − |s (n t ) ) where f 1 and f 2 depend on n t and Γ. An example of the dependence on n t is shown in Fig. 10 for Γ = γ 3 γ 5 for the CG and the BiCGstab2 solvers. f 1 is roughly independent of n t and f 2 decreases rapidly with increasing n t . This behaviour is expected since after a few iterations the first term of Eq. (13) contains most of the signal (and its error) while the second term (and its error) approaches zero. The results were generated on a single configuration at κ loop = 0.166 using 300 stochastic sources. Convergence is achieved after n conv ≈ 480 CG iterations.  Table 7: Gains obtained for various n t , setting N 1 /N 2 = f 1 / f 2 , at κ loop = 0.166, using the BiCGstab2 solver.
The approximate cost is given by, cost ≈ N 1 n t + N 2 n conv .

(A.2)
When TSM is combined with HPE and/or TEA, there are corrections to this formula which the reader can easily work out. Using Lagrange multipliers, we minimize the variance Eq. (A.1) at fixed cost, assuming f 1 to be approximately independent of n t . This yields the optimal values, where f ′ 2 = df 2 /dn t . The right hand sides of Eqs. (A.3) and (A.4) can be computed as functions of n t on one configuration, using a single set of stochastic sources. By finding the intersection between the curve f 1 f 2 /(n conv f ′ 2 2 ) and n t , one can extract n opt t and subsequently determine N 1 /N 2 for this n t . Fig. 11 illustrates this procedure for Γ = γ 3 γ 5 for the CG solver. We read off the values n opt t ≈ 18 and N 1 /N 2 ≈ 30 from this figure, see Table 1.
For our observables we find that when using the optimized TSM parameter values, the two variances within Eq. (A.1) are of similar sizes, f 1 /N 1 ≈ f 2 /N 2 . It also turns out that the gain factor does not critically depend on N 1 /N 2 . For instance, increasing this ratio from the optimal value of 30 found for Γ = γ 3 γ 5 at n t = 18 to the equal cost value N 1 /N 2 = f 1 / f 2 ≈ 35 will increase the final variance by just 3 %, which is statistically insignificant. This suggests an alternative criterium for fixing the parameters: scanning through n t , keeping N 1 /N 2 = f 1 / f 2 fixed, to determine the n t value with the smallest combined variance. We find that following this strategy reduces no gain factor by more than 2 %, relative to the gain achieved using the optimal values, at κ = 0.166, using the CG solver.
It can be seen from Fig. 10 that the convergence is no longer a smooth function of n t when using the BiCGstab2 solver, so that f ′ 2 cannot be determined. The situation is even worse for Γ = ½. In the BiCGstab2 case we have to resort to the alternative method discussed above. Table 7 demonstrates that, using this approach, one can indeed find values n t and N 1 /N 2 for the BiCGstab2 solver that give reasonable gains. The cost was kept fixed to correspond to 300 BiCGstab2 solves to convergence. In this case n conv ≈ 156 is by a factor three smaller than for the CG, however, each iteration is about twice as expensive. The best gain factors are 5.0 for Γ = ½ and 7.3 for Γ = γ 3 γ 5 while for CG, using this method, we are able to achieve very similar factors of 4.9 and 7.9, respectively. The CG algorithm is more robust than BiCGstab2 and gives nearly optimal gains over a wider range of n t values. Note for instance the somewhat erratic behaviour in Table 7 of BiCGstab2 at n t = 10. Moreover, the optimal N 1 /N 2 ratios come out rather large, due to tiny f 2 values, which turns BiCGstab2 less optimal when it is combined with HPE or TEA. Therefore, in the context of TSM, the CG solver is our preferred choice. However, it is possible to combine both solvers, using CG for the truncated solves and BiCGstab2 for running to convergence more efficiently.