Better estimation of real-time correlations

[kbarros]

This greatly reduces artifacts. Is there a reason not to?

Repurposing this issue for a discussion of possible improvements to the estimation of dynamical correlations. We currently have process_trajectory=:symmetrize, analyzed below, and @Lazersmoke is working on a different approach that avoids "periodic wrapping" of correlations on a finite time trajectory.

For the latter, here is a note that builds on @Lazersmoke's ideas:
causal_correlations.pdf

[ddahlbom]

I wonder if the default should be windowing. That seemed to have the greatest effect in the results Sam was reporting at the end of last year. Xiaojian the other day was pointing out the symmetrization doesn't work as well as he'd like in some cases.

[kbarros]

Could be. Shall we collect data for the various S(q,w) results for some simple, prototypical model? How about we compare LSWT, small T classical dynamics with/without :symmetrize in current Sunny, and with @Lazersmoke's upcoming windowing scheme?

[Lazersmoke]

Two comments on this. Firstly, is the implementation of this correct? It seems to me like this for loop is only correct up to t/2, since it is reading and writing from samplebuf on every iteration:

Sunny.jl/src/SampledCorrelations/CorrelationSampling.jl

Line 67 in 3af28c5

for t in 1:nsteps

Secondly, the effect of this option is exactly the same as using the correlations of the time-reversed trajectory (of each sampled trajectory) as an additional estimator of the true correlations. This is certainly not correct in the general case (see the test_correlation_sampling test case on the asymmetric-rebased branch).

There is some reason to suspect that the correlations may be time-symmetric when the system is in equilibrium. So having :symmetrize as an option to make Sunny output correlations which are exactly time-symmetric is useful! But I would not make it default because it's certainly unphysical for some out of equilibrium systems.

[kbarros]

Correctness of statistics from time-reversed trajectories seems contingent on:

1. The assumption that spin configurations are sampled from a good thermal equilibrium.
2. The time-reversal symmetry of the dynamics.

I think 1. is a valid assumption for the large majority of use cases (although @Lazersmoke's research is a notable exception). I think I'd be OK with assuming equilibrium by default (with an "opt out" option) if that were to provide advantages in the common case. I think 2. could be more concerning -- if the magnetic ordering breaks chiral symmetry, is it possible that the dynamics is not time-reversal symmetric? Perhaps @hlane33 also has insight.

[ddahlbom]

To @Lazersmoke's first point, I think they are right that the implementation is wrong. The combination of forward and reverse trajectories is incorrectly weighted after passing the midpoint. If we want to avoid having an allocation or a buffer, that loop can just be run to the midpoint and then the second have can be filled with the time-reversal of the completed first half (or do the setting of mirrored points simultaneously).

To @kbarros's point, we did discuss a simple model some time last year, I think in November. I believe our main conclusion, with respect to just eliminating artifacts due to the discontinuity at the periodic extension, was that a simple cosine window was preferable to symmetrization. That said, it looks like the implementation of symmetrization was incorrect, so it would be worth revisiting.

[kbarros]

After playing with the math, I agree now that :symmetrize seems dubious.

[kbarros]

I spent some more time on this, and I think there could be something interesting here. See attached note. After averaging over an arbitrary phase, a factor of $\sqrt{2}$ appears, and it seems possible to get the exactly correct structure factor?
symmetrized_signal.pdf

[Lazersmoke]

As can be seen by the sqrt(2), the correlations depend non-linearly on the trajectories. So the symmetrized procedure actually involves a linear combination of the four correlations <AB> <A'B'> <A'B> and <AB'> where X' is the time-reversed trajectory of X. So while in the simple cosine case, things are time-reversal symmetric enough (at least after phase averaging) to make this combination coincide with the right answer, this doesn't extend to the realistic cases (can't extend by linearity because it's not linear)

[kbarros]

For posterity, here is an extension of the note for calculating the cross correlation between two different real signals A(t) and B(t). It appears that the technique will still work in this more general case, and an overall correction factor of sqrt(2) is still needed.

I suspect that Sam's more recent approach (zero padding to mask unphysical correlations and then cosine windowing to dampen statistical noise) will be overall better, but it would be good to get a detailed comparison before selecting a new default.

symmetrized_signal2.pdf

[Lazersmoke]

In this note, it's shown that, after equilibrium averaging, the quantity c(omega) + c(-omega) can be estimated using the :symmetrize flag. In other words, the time-symmetrized correlations can be computed. (They can also be computed "in post" by symmetrizing the possibly asymmetric correlations, though)

However, this is not the same thing as the real part of c(omega) (due to the spatial Fourier transform). Further, <SxSy> is a perfectly legitimate quantity to consider, mainly for comparison to the quantum case.

The upshot is that this note means that we can document the option as providing "time-symmetrized correlations (on equilibrium average)" :)

[kbarros]

We discussed on Slack. The symmetrized_signal2.pdf note seems mathematically correct, but there are caveats in the physical interpretation:

1. The note assumes that we are talking about time correlations of real signals $A(t)$ and $B(t)$. If one starts with a complex signal $A(t)$, e.g., as obtained from $\langle S_x S_y \rangle = \mathbf{Z}^\dagger S_x S_y \mathbf{Z}$, then it is not true that $A_\omega^* = A_{-\omega}$, which is a property used at the beginning of the note. It's not clear whether this is a problematic limitation in practice, because observable operators must always be Hermitian, leading to real expectation values. For example, $A(t)$ might represent the time signal for $\langle S_x S_y + S_y S_x \rangle$, and this is fine. [Update: I think this concern gets resolved if "symmetrization" means time-symmetrize and complex-conjugate, simultaneously.]
2. Even if one starts from real signals $A(t)$ and $B(t)$, their cross correlation $C(t) = A \star B$ need not be symmetric in time. Symmetrizing $A(t)$ and $B(t)$ leads to symmetric $\tilde{C}(t) \equiv [C(t) + C(-t)]/2$. In Fourier space, this means $\tilde{C}(\omega)$ is purely real, and we have lost all information about the imaginary part of $C_\omega$, which may be important. It is only in the special case that $A(t) = B(t)$ that we do not lose information.

[Lazersmoke]

It turns out to be very important that the window function applied to the time correlations (or the "default" box window) is exactly symmetric. Otherwise you can get imaginary valued or negative intensities. I've implemented this here:

ae3277b (lives in #217 )

The subtlety is that xs[1] is the equal-time correlations, but xs[end] and xs[2] are un-equal time correlations, so you can't use a completely naive "cosine squared window"; it has to be shifted by one sample!

[Lazersmoke]

Just another data point I found: QuantumOptics.jl also has the capability to compute correlations and their fourier transform (for certain quantum systems).

Their approach involves computing a trajectory on the interval [0,T], and applying what they call "Wiener-Khinchin theorem", which is to take twice the real part of the 'unilateral fourier transform': $2Re\int_0^T d\tau\ e^{i\omega \tau} \langle A^\dagger(t+\tau) A(t)\rangle$. The idea is that they are only interested in autocorrelations so (at least when $A$ is a normal operator, not sure about otherwise), the $A^\dagger A$ operator product is always real (even in the quantum case). And (apparently) because in their context, the angle brackets denote a $tr(\rho \cdots)$ type average over a time-invariant (steady state by assumption) $\rho$, they further expect exactly time-symmetric correlations.

Because of all of those reasons listed above, the $2Re$ transformation appears to be valid in their use case. Unfortunately, it looks like they don't zero-pad their FFT :(. I took their test case and increased the time resolution to tspan = [0.:0.01:15;], and recomputed both real.(exp_values) (the directly computed direct-time correlations; blue) and real.(ifft(ifftshift(SFFT))) * 50, (the inverse fourier transform of the $2Re$ computed spectrum) with the factor 50 needed to match the equal-time correlation, and plotted them:

As you can see, they don't match. The point is to say that this is subtle; we zero-pad so we won't have this bug, but generally speaking the bugs are difficult to detect.

[kbarros]

This was fixed in #246.