getting smearing to work

[kostrzewa]

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).

[kostrzewa]

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.

[kostrzewa]

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?

[deuzeman]

That is a very useful comparison to have, nice work! What sort of configuration are these force comparisons made on, though?

[kostrzewa]

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?

[kostrzewa]

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.

[kostrzewa]

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...

[kostrzewa]

    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]);

[kostrzewa]

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 :)

[kostrzewa]

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].

[kostrzewa]

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.

[kostrzewa]

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.

[urbach]

the force of the smeared part should be independent of the monomials at all, its just the chain rule, isn't it?

[kostrzewa]

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.

[kostrzewa]

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

[kostrzewa]

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

[kostrzewa]

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.

[deuzeman]

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.

[kostrzewa]

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...)

[kostrzewa]

[...] 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.

[deuzeman]

Just to be concrete, you are referring to Sigma' rather than Sigma, correct?

Yes, I was.

[urbach]

Hey guys, merry christmas!

[kostrzewa]

:) Merry Christmas Carsten and Albert :)

[deuzeman]

:D A merry Christmas to both of you as well! :christmas_tree:

[deuzeman]

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...

[kostrzewa]

I give up... I've checked the calculation and I think I now match P&M everywhere so...

[kostrzewa]

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?

[urbach]

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...!?

[kostrzewa]

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.

[urbach]

are we using the correct basis for the generators?

[kostrzewa]

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).

[deuzeman]

A happy New Year to you as well! I'm quite happy that you did the explicit BMW vs PM check, because that was on my to-do list. I looked at it before and had concluded that they ended up giving the same result in spite of some superficial differences in the formulation. But then there's always this nagging feeling that you did a shoddy job. Seeing agreement here does exclude a number of possible errors, even though I guess there's still much shared code.

For now, I will try to confirm my earlier conclusion about the deviation between PM and BMW for my Abelian test case. It would seem I was sloppy there. If they do agree with each other, there may still be the mathematical disagreement. In that case, maybe there's an issue with our interpretation of the formulae?

[kostrzewa]

Regarding your worries about Z, according to 3.32 when it is eventually used the extra factor of U^dagger is cancelled off by the multiplication from the left by U.

From the same equation I can see that the only difference in the final force is a multiplication from the left by U in the first term, by U^dagger from the right in the second term and I guess by U in the third and final term.

[deuzeman]

The P&M paper should make sense one way or the other, or I'm sure Istvan would have mentioned something. But we may be interpreting it wrongly. And then I am also thinking of something like a subtle issue in the translation between different representations, such that what we use as Sigma is not what appears in their derivation. Not sure where that would be happening, though. We're mainly recycling routines from the rest of the code here, so it's hard to imagine we're mixing up conventions or something.

I don't think I've pushed all the testing modifications to Github, mainly because it produces a pretty mangled executable. But I can make the branch available.

The issue that I saw should not be related to the projection, though. Since I was working with diagonal matrices everywhere, the projection simply amounts to taking the imaginary component of the matrix.

[kostrzewa]

And then I am also thinking of something like a subtle issue in the translation between different representations, such that what we use as Sigma is not what appears in their derivation.

I also doubt this should be a problem because our projections are inverses of eachother, so it doesn't really matter which particular form of the generators we use...

[deuzeman]

Regarding your worries about Z, according to 3.32 when it is eventually used the extra factor of U^dagger is cancelled off by the multiplication from the left by U.

True, that factor does get cancelled. But I'm not sure the factor V^dagger in Sigma (3.24) can be cancelled in a similar manner when Sigma appears rather non-trivially in the definition of Z. Mind you, I did at some point convince myself that there was no difference between the two definitions. So I was either overlooking something then, or I'm doing it now... ^_^

From the same equation I can see that the only difference in the final force is a multiplication from the left by U in the first term, by U^dagger from the right in the second term and I guess by U in the third and final term.

I think the sandwiching works out in the end, it should be a convenient rewriting of P&M's expression anywhere, except possibly for Sigma in Z and, as far as I'm aware, for the two terms with a negative prefactor in the staple derivative. For the latter, it seems there's some hermitian conjugate being taken and I think the difference between the two definitions is cancelled in the projection.

[kostrzewa]

But I'm not sure the factor V^dagger in Sigma (3.24) can be cancelled in a similar manner when Sigma appears rather non-trivially in the definition of Z

Yes, I agree... also, the post-multiplication of sigma by V in 3.32/3.33 is essential to recover exp(A) [through post-multiplication with U^dagger]. This V should therefore better not be involved in any cancellations!

[urbach]

Hey Albert, isn't it possible to backtrace where the non-zero eights element came from?

[kostrzewa]

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.

Which definition of lambda_3 are you using? When I run such a matrix to exposu3_inplace this results in the unit matrix and when I run it though cayley hamilton I get a diagonal matrix proportional to the unit matrix. I might be using the wrong lambda though...

~ lambda_3
[ + 1.550000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 1.550000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]

exp_in_place(~lambda_3) [projects onto gell-mann and then exponentiates that]
[ + 1.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, + 1.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, + 1.000000000000 - 0.000000000000 * i]

~cayley hamilton(~lambda_3)
[ + 0.020794827803 + 0.999783764189 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, + 0.020794827803 - 0.999783764189 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, + 1.000000000000 - 0.000000000000 * i]

[kostrzewa]

Which definition of lambda_3 are you using?

I guess you used the traceless antihermitian version, correct?

[kostrzewa]

Doing the same thing I first of all see that Cayley hamilton breaks down for this... But exposu3_inplace does what you say. However, I see no such deviation in the eighth component (with my P&M code):

P&M

~ lambda_3 (TA)
[ - 0.000000000000 + 1.550000000000 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 1.550000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]

 gauge field at [1][1]
[ + 0.020794827750 + 0.999783764185 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, + 0.020794827750 - 0.999783764185 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, + 1.000000000000 - 0.000000000000 * i]

force:
[DEBUG] Comparison of force calculation at [1][1]!
         numerical force <-> analytical force <-> analytical force + smeared 
    [0]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [1]  +0.000000000000 <-> +0.000000000000 <-> -0.000000000000
    [2]  -4.873255619486 <-> -5.236968500253 <-> +4.080262758119
    [3]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [4]  +0.000000000000 <-> +0.000000000000 <-> -0.000000000000
    [5]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [6]  +0.000000000000 <-> +0.000000000000 <-> -0.000000000000
    [7]  +0.000000000000 <-> +0.000000000000 <-> -0.000000000000

[kostrzewa]

I first of all see that Cayley hamilton breaks down for this

this is, of course, for the Cayley Hamilton as used in P&M where it returns exp(iM)...

[kostrzewa]

Ooops... I made a big mistake.. sorry! I had the P&M commit in what I though was my BMW branch.. I only noticed after the last message.. I can confirm that the BMW version of the code produces the weirdness observed by Albert:

BMW

[ - 0.000000000000 + 1.550000000000 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 1.550000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]

[ + 0.020794827750 + 0.999783764185 * i, - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, + 0.020794827750 - 0.999783764185 * i, - 0.000000000000 - 0.000000000000 * i]
[ - 0.000000000000 - 0.000000000000 * i, - 0.000000000000 - 0.000000000000 * i, + 1.000000000000 - 0.000000000000 * i]
# Initialised monomial of type GAUGE, no_monomials= 1
# monomial id 0 type = 2 timescale 0
# Computed plaquette value: 0.959199784490.
#
# Starting trajectory no 1
called gauge_heatbath for id 0 1 energy0 = 488.040850
# Time for numerical GAUGE monomial derivative: 7.811116e-02 s
# Time for analytical GAUGE monomial derivative: 2.783700e-05 s
[DEBUG] Comparison of force calculation at [1][1]!
         numerical force <-> analytical force <-> analytical force + smeared 
    [0]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [1]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [2]  -4.873255619486 <-> -5.236968500253 <-> -5.324035854924
    [3]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [4]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [5]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [6]  +0.000000000000 <-> +0.000000000000 <-> +0.000000000000
    [7]  +0.000000000000 <-> +0.000000000000 <-> -1.147610951801

[kostrzewa]

And hence there is also no agreement between P&M and BMW:

BMW

    [0]  -0.262370019755 <-> +0.107905120571 <-> +0.129676042557
    [1]  -0.887668051064 <-> +0.467106242213 <-> +0.471951073001
    [2]  +0.241326313244 <-> -0.160306431233 <-> -0.182220618397
    [3]  -0.225912856422 <-> +0.246215513575 <-> +0.252953544621
    [4]  +0.404052889280 <-> -0.250874926745 <-> -0.240528280807
    [5]  +0.635294225049 <-> -0.407611167750 <-> -0.408778522129
    [6]  -0.295563387454 <-> +0.302704762536 <-> +0.288898460358
    [7]  +0.847369193480 <-> -0.603752630725 <-> -0.591925958260

P&M

    [0]  -0.262370195969 <-> +0.107905141341 <-> -0.110081222496
    [1]  -0.887668170435 <-> +0.467106247189 <-> -0.127970716785
    [2]  +0.241326529249 <-> -0.160306450699 <-> -0.089983558874
    [3]  -0.225912907581 <-> +0.246215547881 <-> -0.045263761964
    [4]  +0.404052951808 <-> -0.250874930380 <-> +0.064943694746
    [5]  +0.635294054518 <-> -0.407611191183 <-> +0.018747520689
    [6]  -0.295563279451 <-> +0.302704711049 <-> +0.006915465669
    [7]  +0.847369142321 <-> -0.603752667446 <-> +0.042348747294

As you can see the smearing is the same but the force computation is not compatible. It's a shame neither of them is correct though at this stage...

My next step was to try to introduce a different projector which uses the standard form for the Gell-Mann matrices and see what happens.

[kostrzewa]

Sorry again about the confusion.. I don't know how that commit snuck in there...

[kostrzewa]

Maybe we should ask Istvan whether he know details about the development of the P&M implementation in 'his' code?

[deuzeman]

It's good that you see the same issue, but indeed a shame that the P&M formalism doesn't produce the right answer either. Still, it does a better job than the rather random mixing of BMW...

First, with regards to your remarks on the Cayleigh-Hamilton method... It wasn't quite clear to me if you think there's an issue there in the end. I have seen it produce correct results, when comparing to Urs' stout smearing code and to Matlab's results for exponentiation. But it's been written to work with BMW's conventions, so it actually calculates exp(A) for A being a traceless antihermitian matrix. Was that accounted for?

Of course, the values of B1 and B2 cannot be checked easily in that manner, so they may still be buggy somehow. They involve some subtle cancellations, however, leading to obviously problematic values when I had a sign wrong earlier. There could, however, be an overall phase missing for those matrices. Maybe something to look into as well?

[kostrzewa]

First, with regards to your remarks on the Cayleigh-Hamilton method... It wasn't quite clear to me if you think there's an issue there in the end. I have seen it produce correct results, when comparing to Urs' stout smearing code and to Matlab's results for exponentiation. But it's been written to work with BMW's conventions, so it actually calculates exp(A) for A being a traceless antihermitian matrix. Was that accounted for?

Yes, I changed it to compute exp(iQ) and I now store Q in place of A and I've taken care to change all factors of (-i)^2 and (-i) that you introduced. This is why it blew up when I exponentiated the matrix further above. (so it doesn't work as a replacement for exposu3_inplace anymore)

Of course, the values of B1 and B2 cannot be checked easily in that manner, so they may still be buggy somehow. They involve some subtle cancellations, however, leading to obviously problematic values when I had a sign wrong earlier. There could, however, be an overall phase missing for those matrices. Maybe something to look into as well?

Well, I checked the calculation line by line w.r.t. to the P&M paper and I can't find fault with it.

I have found the stout implementation in Istvan's code and am currently trying to decode what's done. Luckily it follows the P&M paper to the letter. (not even optimizing multiplications) (/homea/hch02/hch02c/ForAbdou/UpdateCode.JUGENE In Juelich)

[kostrzewa]

Yes, I changed it to compute exp(iQ) and I now store Q in place of A and I've taken care to change all factors of (-i)^2 and (-i) that you introduced. This is why it blew up when I exponentiated the matrix further above. (so it doesn't work as a replacement for exposu3_inplace anymore)

oh wait, but then it computes the wrong B1 and B2, no? I mean, I changed exponent_from_coefficients.. Oh wait, no it computes it correctly since I store Q!

[deuzeman]

Potentially, yes. But I think they may be derivatives of IQ, rather than Q. I think I looked at it and concluded it was fine. But since something must be wrong somewhere... On 4 Jan 2013 12:03, "Bartosz Kostrzewa" notifications@github.com wrote:

Yes, I changed it to compute exp(iQ) and I now store Q in place of A and I've taken care to change all factors of (-i)^2 and (-i) that you introduced. This is why it blew up when I exponentiated the matrix further above. (so it doesn't work as a replacement for exposu3_inplace anymore)

oh wait, but then it computes the wrong B1 and B2, no?

— Reply to this email directly or view it on GitHubhttps://github.com/etmc/tmLQCD/issues/198#issuecomment-11879731.

[urbach]

Hey! I think you are making really good progress! Thanks for all of this! Not having looked into the details of the code, but where does the non-zero eights element show up exactly. I would think that this can be traced back? And then its maybe easier to understand whats going wrong (which was my original intention by proposing this kind of check).

I'll try to implement a rational monomial while you debug smearing, unless I can be of help somewhere...!