etmc
getting smearing to work
I've invested a bit of time into helping Albert in getting this rather complex but beautifully written piece of code working properly. I would like to use this issue to somewhat track my progress.
I started by turning the numerical derivative into a symmetric one to remove linear discretization errors. I then added a real numerical derivative for the Wilson plaquette gauge action to the gauge monomial. One trajectory resulting in a 100% change in the plaquette measurement only has a difference in dH of 1e-8 or so and acceptance is the same when compared to the trajectory computed with the analytical derivative.
Adding support for smearing was a bit trickier so I add a stout smearing control with hardcoded parameters in the gauge monomial derivative itself and destroy it afterwards. I also explicitly disabled the computation and accumulation of smearing force terms in update_momenta and stout_stout_smear.
With these constructs in place and the same parameters in the input file I am now able (apparently) to perform pure-gauge smeared HMC with 100% acceptance. (there is an input file in the branch https://github.com/kostrzewa/tmLQCD/tree/numerical_deriv_num_gauge_monomial ) The code still fails to produce acceptance for more than one iteration of stout smearing which is a bit strange.
Maybe the buffer swapping doesn't work properly? Or maybe there needs to be a separation between the updating of the gauge energy and the various acceptance and heatbath functions? Right now it is unfortunately necessary to set g_update_gauge_energy explicitly many times and hope that measure_gauge_action is called with the correct gauge field next time around... This really should be changed but since we'll never calculate the gauge energy on a smeared gauge field, this can wait as long as we make sure it's always called on g_gf.
Finally a comment on the numerical derivative. I think the type of thing implemented by Albert (rotations on one gauge link only, 8 components) is a great sanity check tool to implement on for a more general per-monomial basis. (call analytical derivative, compute, check for a flag, compute numerical derivative for this one gauge link, display output, at debug level 5, say) I'm not sure how it would be best integrated into the code though to be practical and inconspicuous. Since it scales with 64*VOLUME^2, it is unrealistic to do a full numerical derivative for anything other than the gauge monomial (and even here only for 2^4 volumes).
I think the type of thing implemented by Albert (rotations on one gauge link only, 8 components) is a great sanity check tool to implement on for a more general per-monomial basis.
This needs a little bit of work to take into account parallelization. In that case it would also be better to compute something at a boundary and something in the bulk as an additional check.
Testing smearing forces
I've now re-enabled the calculation of the force contribution of the smearing but I don't accumulate it in df at the end. I call the monomial derivative with df0 set to df (no remapping of df0!) and then explicitly call smear_forces with df as an input.
smearing_control_monomial[s_type]->force_result[1][1] then contains the following compared to df in the first trajectory, the first time the forces are computed:
# Time for GAUGE monomial derivative: 5.819154e-02 s
[DEBUG] Comparison of force calculation at [1][1]!
Involves monomials: GAUGE
Numerical total force <-> smear forces
[0] +0.677441909147 <-> +0.674610818420
[1] +0.651245056815 <-> +0.622717380971
[2] -0.650447663020 <-> -0.645114865239
[3] -0.384530363817 <-> -0.366998593273
[4] +0.027326740337 <-> +0.057354221129
[5] -0.583583420166 <-> -0.604351922979
[6] +0.865918178761 <-> +0.860550664196
[7] -0.944989812979 <-> -0.944308827867
Soon, however (still in the same trajectory), the two have diverged subtantially:
[DEBUG] Comparison of force calculation at [1][1]!
Involves monomials: GAUGE
Numerical total force <-> smear forces
[0] +2.188621795085 <-> +2.513886660328
[1] +2.763872868172 <-> +2.223844539664
[2] -2.256187411831 <-> -2.343240983808
[3] -0.849516123935 <-> -0.814242697326
[4] +0.072232415960 <-> +1.132103245325
[5] -2.549349795800 <-> -2.980789873106
[6] +3.128423037424 <-> +3.023725888304
[7] -4.163150993008 <-> -4.485607838646
If I can trust what I did then I would suspect that there is some buffer or something else which is not properly cleared between calls to smear_forces, therefore accumulating an error.
How large is the contribution to the force coming from the smearing expected to be?
That is a very useful comparison to have, nice work! What sort of configuration are these force comparisons made on, though?
That is a very useful comparison to have, nice work! What sort of configuration are these force comparisons made on, though?
I guess I don't really understand the question. Since I have only one monomial and this one is smeared there is one relevant smearing and one smeared gauge field in update_momenta. This smeared gauge field should be, by construction, the same one that is used in the monomial derivative (as there I memmove from g_gf, rotate, then smear with a local control with the same parameters, so if I smeared g_gf right there and then, I would obtain the same smeared gauge_field as here in update_momenta)
for (int s_ctr = 0; s_ctr < no_relevant_smearings; ++s_ctr)
{
int s_type = relevant_smearings[s_ctr];
smear(smearing_control_monomial[s_type], g_gf);
ohnohack_remap_g_gauge_field(smearing_control_monomial[s_type]->result);
g_update_gauge_energy=1;
g_update_rectangle_energy=1;
for(int k = 0; k < no; k++)
{
if (monomial_list[ mnllist[k] ].smearing == s_type)
{
if(monomial_list[ mnllist[k] ].derivativefunction != NULL)
{
#ifdef MPI
atime = MPI_Wtime();
#else
atime = (double)clock()/(double)(CLOCKS_PER_SEC);
#endif
// as df0 has not been remapped, this writes directly into df
monomial_list[ mnllist[k] ].derivativefunction(mnllist[k], hf);
#ifdef MPI
etime = MPI_Wtime();
#else
etime = (double)clock()/(double)(CLOCKS_PER_SEC);
#endif
}
}
}
smear_forces(smearing_control_monomial[s_type], df);
[ output comparison ]
}
ohnohack_remap_g_gauge_field(g_gf);
That should be legitimate, no?
Okay, found one bug just need to figure out how to fix it. In stout_smear_plain, when multiple iterations are used the result blows up into NaN. I'm presuming, again, a fault in the buffer swap or so. If I explicitly construct the control with "compute force terms" enabled (and then simply discard the intermediate stuff from the call to stout_smear_with_force_terms), multiple iterations work just fine and I have acceptance.
I looked into your concern regarding "which configuration the force is computed on" and moved the comparison into the gauge monomial. Here, for one particular component, direction and gauge link I rotate, smear and then smear_forces. The result does not agree with the numerical one. I think that should be a good sign that the force terms are wrong, no? Although I think what I did before was more correct...
stout_control* control = construct_stout_control(1,2,0.18);
for(int x = 0; x < VOLUME; ++x)
{
for(int mu = 0; mu < 4; ++mu)
{
xm=(double*)&hf->derivative[x][mu];
control->smearing_performed = 0;
for (int component = 0; component < 8; ++component)
{
double h_rotated[2] = {0.0,0.0};
printf("Rotating at %d %d in component %d\n",x,mu,component);
for(int direction = 0; direction < 2; ++direction)
{
link=&rotated[direction][x][mu];
memmove(&old_value, link, sizeof(su3));
/* Introduce a rotation along one of the components */
memset(ar_rotation, 0, sizeof(su3adj));
ar_rotation[component] = epsilon[direction];
exposu3(&mat_rotation, &rotation);
_su3_times_su3(rotated[direction][x][mu], mat_rotation, old_value);
stout_smear(control, rotated[direction]);
if(x == 1 && mu == 1 && direction == 1 && component == 1 ) {
stout_smear_forces(control,df);
fprintf(stderr, "[DEBUG] Comparison of force calculation at [1][1]!\n");
fprintf(stderr, " smear forces <-> numerical total force\n");
fprintf(stderr, " [%d] %+14.12f <-> ", component, control->force_result[1][1].d1); //*/
}
// compute gauge action
g_update_gauge_energy = 1;
g_update_rectangle_energy = 1;
h_rotated[direction] = -g_beta*(mnl->c0 * measure_gauge_action(control->result)); //rotated[direction]));
// reset modified part of gauge field
memmove(link,&old_value, sizeof(su3));
} // direction
// calculate force contribution from gauge field due to rotation
xm[component] += (h_rotated[1]-h_rotated[0])/(2*eps);
if( x == 1 && mu == 1 && component == 1 )
fprintf(stderr, "%+14.12f\n", df[1][1].d1); //*/
} // component
} // mu
} // x
free_stout_control(control);
return_gauge_field(&rotated[0]);
return_gauge_field(&rotated[1]);
Just as a side note I can now confirm that I have acceptance even with 4 iterations of stout smearing. For this small volume this pushes the plaquette up to 0.98 though :)
Although I think what I did before was more correct...
Indeed, the very first time the force is computed the result from smear_forces in gauge_monomial.c is further from the numerical result. In my opinion this makes perfect sense because here the forces are computed on a rotated and then smeared gauge_field and not on one that has just been smeared.
[DEBUG] Comparison of force calculation at [1][1]!
smear forces <-> numerical total force
[1] +0.949474738841 <-> +1.300266217186
[DEBUG] Comparison of force calculation at [1][1]!
Involves monomials: GAUGE
Numerical total force <-> smear forces
[0] +1.300266217186 <-> +1.376853348262
[1] +1.643219286507 <-> -0.445783308699
[2] +5.559405138911 <-> +5.240163473078
[3] -5.011576001834 <-> -4.214708159595
[4] -3.231603625409 <-> -1.894558971907
[5] -0.102085107301 <-> -1.403325081582
[6] +3.653467098275 <-> +0.851817872119
[7] -7.133829001305 <-> -7.357334898408
In any case I think we need to investigate the force computation now which will inevitably involve going through the derivation...
UPDATE
By the way, the second output still comes from update_momenta. I didn't actually notice the massive deviations in [1], [4], [5] and [6].
In any case I think we need to investigate the force computation now which will inevitably involve going through the derivation...
Well, first I will go through the code line-by-line to check if I don't notice something that has gone by.
One thing that perplexes me a little is the lack of any mention of beta in the computation of the force terms. The gauge field is stored independently of beta which is why in the computation of the gauge derivative there is an explicit multiplication by beta (in the form of "factor"). Does rho implicitly contain beta? Anyway, I'll read the papers and notes and then I'm sure I'll understand.
the force of the smeared part should be independent of the monomials at all, its just the chain rule, isn't it?
the force of the smeared part should be independent of the monomials at all, its just the chain rule, isn't it?
The output here
[DEBUG] Comparison of force calculation at [1][1]!
smear forces <-> numerical total force
[1] +0.949474738841 <-> +1.300266217186
is wrong because it is computed on a rotated gaugefield and cannot be correct. Further, the comparison I'm making in update_momenta is also wrong because the derivative that is fed into smear_forces already contains the smearing term implicitly. I will compute the derivative of the gauge monomial seperately analytically and call smear forced on this temporary derivative to see if any differences are present.
I changed the code a bit more but I still have my doubts... I first compute the numerical derivative as explained before (the gauge field is the smeared one). Then I compute the analytical derivative but I remap df0 first, the gauge field is still mapped to the smeared field. Now I call smear_forces on this derivative.
Does the numerical derivative really capture the complete influence of the smearing or is there somethig missing, something too much? Since I get acceptance I would think that everything is correct,
# Time for GAUGE monomial derivative: 3.404317e-01 s
# Time for GAUGE monomial derivative: 8.992700e-05 s
[DEBUG] Comparison of force calculation at [1][1]!
numerical force <-> analytical force <-> analytical force + smeared
[0] -1.399133145696 <-> -1.614094684864 <-> -1.581404238899
[1] +1.379116918088 <-> +1.273875681870 <-> 1.230786767312
[2] +3.851927692722 <-> +3.154854008718 <-> 3.174869727485
[3] -2.446024618052 <-> -1.753064895060 <-> -1.745692061800
[4] -0.683806388224 <-> -0.141420359739 <-> -0.104222017625
[5] +2.037271710265 <-> +2.023458429842 <-> 1.976318650656
[6] -0.187620292991 <-> -0.757265165433 <-> -0.818135956147
[7] -2.591452161482 <-> -1.657051994768 <-> -1.692047414142
# Time gauge update: 1.118730e-04 s
Reducing rho subtanitally seems to indicate that the numerical derivative is wrong after all... actually no.. forget what I just said!! Just had a parameter not upated. For rho 0.0001:
[DEBUG] Comparison of force calculation at [1][1]!
numerical force <-> analytical force <-> analytical force + smeared
[0] -1.716202169177 <-> -1.716295097285 <-> -1.716294988623
[1] +1.106646996618 <-> +1.106390433190 <-> +1.106390261847
[2] +2.880228134927 <-> +2.879208553427 <-> +2.879208658993
[3] -1.236969762886 <-> -1.235798700022 <-> -1.235798717189
[4] +0.409221449615 <-> +0.410269441569 <-> +0.410269551942
[5] +2.118844865606 <-> +2.118757717394 <-> +2.118757550812
[6] -1.362196374544 <-> -1.363141425826 <-> -1.363141647349
[7] -1.145151750848 <-> -1.143913614477 <-> -1.143913783982
# Time gauge update: 1.133960e-04 s
Without smearing (rho=0.0) everything agrees up to discretisation errors:
# Time for GAUGE monomial derivative: 1.753207e-01 s
# Time for GAUGE monomial derivative: 8.376600e-05 s
[DEBUG] Comparison of force calculation at [1][1]!
numerical force <-> analytical force <-> analytical force + smeared
[0] -1.206807453968 <-> -1.206807453026 <-> -1.206807453026
[1] +1.138477345819 <-> +1.138477345366 <-> +1.138477345366
[2] +2.416553850537 <-> +2.416553851685 <-> +2.416553851685
[3] -1.283008155895 <-> -1.283008157294 <-> -1.283008157294
[4] -0.301118627988 <-> -0.301118628069 <-> -0.301118628069
[5] +2.658250213017 <-> +2.658250212384 <-> +2.658250212384
[6] -1.354716569324 <-> -1.354716569960 <-> -1.354716569960
[7] -1.729094647374 <-> -1.729094647107 <-> -1.729094647107
the force of the smeared part should be independent of the monomials at all, its just the chain rule, isn't it?
Hi, on the train ride to Luxembourg I read and tried to fully understand the Peardon and Morningstar paper. As far as I did it is not true that the smearing contribution to the force will be independent of the monomial(s) that use(s) smearing. Every part of the force computation implicitly contains Sigma' (through lambda and gamma) and therefore also contains details of the monomial(s).
The first iteration in the force computation is the derivative of the terms of the action that include fat links with respect to the fattest link. This is exactly the same as the standard derivative computation just evaluated on the fattest gauge configuration. The further iterations (there are at least two in total for one iteration of stout smearing) then use this force and the smearing (which types of staples are used in the fattening) definition to compute the full force iteratively. (equation 75 in the P&M paper). Therefore in the case of a pure gauge smeared action like the one the numerical derivative hmc computed, beta enters in the first step implicitly through what you (or Hasenbusch?) called the force factor in the code.
As far as I understand for a fully smeared pure gauge action, equation 41 in P&M would not contain a thin link contribution at all. (It would only enter at the lowest level (\sigma_0) of the analytical force computation, see eq 75 where the LHS would be \sigma_0 and U in the RHS would correspond to thin links.)
Note also that therefore I am now confident that the numerical derivative does the right thing with respect to smearing:
Using V for the smeared gauge field: I compute dS/dU (eq 40) which for this action corresponds to dS_st/dV dV/dU where dV/dU would given by "( V(U+dU)-V(U-dU) )/(2eps)" (formal considerations required to carry out the chain rule "multiplication" )
or alternatively
(S[ V(U+dU) ] - S[ V(U-dU) ]) / 2*eps
this computes the full derivative because dS/dV would be (S[V+dV] - S[V-dV]) / 2*eps which is not the same! In that case only one component of the smeared field would be rotated while in the above definitions the rotation in U propagates to neighbouring links as it corresponds to the steps: rotage gauge field, smear, compute derivative in that order.
For an action in which only the fermionic part is smeared, there would of course still be the standard thin link gauge derivative and no fat link gauge derivative. The fat link derivative would come in through the fermion terms only.
I believe I now have a good enough understanding to attempt to follow the force computation in the code and spot problems. I will just have to read the BMW and Albert's notes carefully when I get back from holidays.
I also think that I can actually compute the iterations numerically by follwing eqn 75. To do this I would first compute the numerical derivative with respect to the smeared field (so smear first, rotate a component in the smeared field and then evaluate the derivative) and then feed this to the final force computation.
Every part of the force computation implicitly contains Sigma' (through lambda and gamma) and therefore also contains details of the monomial(s).
This is true, but the point is that all of the information has already been encoded within the single input force vector. The modifications to that vector made by the smearing routine are completely independent of the details of its production. For that reason, there should really be no entanglement between the code for the specific monomials and the smearing of the forces. It's not clear to me why you implemented the force calculation the way you did.
Note also that therefore I am now confident that the numerical derivative does the right thing with respect to smearing:
Using V for the smeared gauge field: I compute dS/dU (eq 40) which for this action corresponds to dS_st/dV dV/dU where dV/dU would given by "( V(U+dU)-V(U-dU) )/(2eps)" (formal considerations required to carry out the chain rule >"multiplication" )
or alternatively
(S[ V(U+dU) ] - S[ V(U-dU) ]) / 2*eps
this computes the full derivative because dS/dV would be (S[V+dV] - S[V-dV]) / 2*eps which is not the same! In that case only one component of the smeared field would be rotated while in the above definitions the rotation in U propagates to neighbouring links as it corresponds to the steps: rotate gauge field, smear, compute derivative in that order.
For an action in which only the fermionic part is smeared, there would of course still be the standard thin link gauge derivative and no fat link gauge derivative. The fat link derivative would come in through the fermion terms only.
Agreed. Its physically more correct to think of the smearing as a part of the definition of the operators appearing in the action, rather than a procedure applied to the gauge fields. What we care about is the change induced in some scalar product of those operators, given a small change in the underlying gauge field. The correct point to act with a small rotation if we want to know the MD forces is therefore always the original gauge field, rather than what amounts to some arbitrary stage half-way the calculation of the operator.
(formal considerations required to carry out the chain rule "multiplication" )
This part is actually much better discussed in the BMW paper. Due to the constraints imposed by the group structure, one has to explicitly define what one means by the multiplication. Unfortunately, because life would have been much simpler with a proper factorization of the different contributions...
I also think that I can actually compute the iterations numerically by following eqn 75. To do this I would first compute the numerical derivative with respect to the smeared field (so smear first, rotate a component in the smeared field and then evaluate the derivative) and then feed this to the final force computation.
While some procedure like this should work to get the different stages out, I'm not sure how much we'd gain. Since equation 75 in P&M and its equivalent 3.35 in BMW are iterative and one level of stout smearing doesn't work yet, we can first try to work that single level out. In that case, the value of Sigma is taken directly from the monomial and is known to be correct. It would be very beneficial if we could find some way to segregate the different contributions to the force calculation beyond that, though.
This is true, but the point is that all of the information has already been encoded within the single input force vector. The modifications to that vector made by the smearing routine are completely independent of the details of its production.
Yes, I understand and agree!
For that reason, there should really be no entanglement between the code for the specific monomials and the smearing of the forces. It's not clear to me why you implemented the force calculation the way you did.
Only out of convenience because I know that I can't do a numerical computation for anything other than the gauge monomial. It seemed easiest to explicitly disable force smearing in update_momenta and compute the entire derivative in the monomial. (rather than adding a complicated set of conditionals to the steps in update_momenta which call the derivatives) And since I need to rotate first and only then smear, this approach would have been quite clumsy in update_momenta, where I would have to pass the smearing control [or at least the parameters] through to the derivative (thereby needing to adjust the derivativefunction interface and so on...)
[...] and one level of stout smearing doesn't work yet, we can first try to work that single level out.
That was the plan but in a sufficiently general way to make it work for multiple iterations too.
In that case, the value of Sigma is taken directly from the monomial and is known to be correct.
Just to be concrete, you are referring to Sigma' rather than Sigma, correct?
It would be very beneficial if we could find some way to segregate the different contributions to the force calculation beyond that, though.
I absolutely agree, of course.
Just to be concrete, you are referring to Sigma' rather than Sigma, correct?
Yes, I was.
I made some progress. While I hope I can push this a little further soon, I'm going to be somewhat occupied in the coming days by personal matters. That's why I figured I should post what I suspect now. Sorry if it's a little rambling, but I'm afraid this is the best I can do for now. If nothing else, it can function as a reminder to myself.
Since we know the numerical derivative is correct even for smeared code, I figured I'd try to explicitly check a very simple case. First make sure the formulae as we understand them produce the numerical result, then make sure the code replicates the formulae. Any errors would have to show up along the way.
My test case was a unit gauge field, but with the [0][0] component set to the exponential of a matrix proportional to lambda_3. That makes for a diagonal SU(3) element, with the last diagonal element equal to unity. I calculated the forces for a pure gauge simulation, with a single level of stout smearing on that gauge monomial.
This setup produces an immediate deviation on the first element of the forces, but in an interesting way. The force calculated for an unsmeared field is only non-zero for the third entry, which makes intuitive sense. The force as calculated by the monomial on the smeared field has the same structure. And so does the force that is numerically calculated for the smeared action. But the analytical calculation of the smeared force has non-zero entries in the third and eighth entry, so it's mixing with the other diagonal generator of SU(3).
After poking around in the code and going through some calculations with Matlab, I'm now suspecting that there may be a problem in the definition of Z as given by BMW. I may be overlooking something or having a mistake somewhere in my calculations, but it doesn't seem to have the right symmetry properties. Note that its definition, except for the added factor of U^\dagger, is identical to that of \Lambda in the P&M paper. But in the latter, \Sigma is defined without an application of the smeared gauge field. If I use that definition and follow P&M, the first part of the force corrections (up to the staple derivatives) produces a result just in the third entry. If I follow BMW's definition, things break down already there...
I give up... I've checked the calculation and I think I now match P&M everywhere so...
Are we confident that the P&M paper is correct? Or more likely, are we sure that the projection of the force from the adjoint to the defining representation is the correct thing to do? Maybe we should transpose after the projection?
Sorry, how does this fit with Alberts latest observations? You would say its not the definition of Z?
Albert, the code you used for your latest test, is it available? It should be possible to find out where the non-zero eights element comes from...!?
I will need to look at BMW carefully but I think the factor of V drops out in the actual definition of the smeared force.
We are using the functions already present in tmLQCD and we correct for the factor of -2 after projecting back to the adjoint (from trace_lambda).