Sampling
Ensemble
s can be generated
using a Generator.
A generator is any object which has a class with a step()
method that takes a configuration as a dictionary and returns a new configuration as a dictionary.
Some generators can measure observables as the step is taken and return them inside the step, in which case the Configurations
need a place to store them, which is created with the inline_observables()
method.
- class supervillain.generator.Generator[source]
For example, a very dumb generator is
class DoNothing(Generator):
def __init__(self):
pass
def step(self, x):
r'''
Just copy the incoming configuration and 'measure' the constant 1
and add it to the configuration.
'''
return x.copy() | {'one': 1}
def inline_observables(self, steps):
r'''
Obviously there is no *actual* reason to compute an observable whose value
is identically 1, but it's useful to show how inline observables work.
'''
return {'one': np.zeros(steps)}
Then, if you
L = supervillain.Lattice2D(5)
S = supervillain.Villain(L, kappa=0.1)
e = supervillain.Ensemble(S).generate(17, DoNothing(), start='cold')
the result will be an ensemble e
that has 17 identical configurations (because DoNothing
makes no updates).
Obviously, for good sampling, DoNothing is the worst imaginable algorithm!
It has an infinitely bad ergodicity problem!
The Villain Formulation
A less dumb algorithm is a local update, which changes only fields in some small area of the lattice.
As an example, we can formulate an update scheme offering localized changes to the \(\phi\) and \(n\) fields in the Villain
formulation.
Picking a site \(x\) at random and proposing a change
for the \(\phi\) on \(x\) and \(n\) on links \(\ell\) which touch \(x\) is ergodic (once swept over the lattice) and satisfies detailed balance so long as we accept the proposal based on the change of action.
The NeighborhoodUpdateSlow
generator implements this update algorithm but suffers from a variety of defects.
First, its implementation makes a lot of calls.
We could ‘unfactor’ the proposal, single-site accept/reject, and sweep of the lattice for some benefits in speed.
Moreover, to compute the change in action it evaluates the action for both configurations and takes the difference.
But, we know that’s silly, most of the \(\phi\)s and \(n\)s haven’t changed at all and those links will contribute to both actions equally, giving zero difference.
We could reduce the arithmetic substantially by computing the difference directly.
Finally, for ease of thinking, each proposal
reckons locations relative the origin and therefore moves all the fields around in order to update ergodically.
All of that data movement adds cost, especially as the lattice gets large.
It was easy to write the NeighborhoodUpdateSlow
but we can do better for production.
- class supervillain.generator.villain.NeighborhoodUpdate(action, interval_phi=3.141592653589793, interval_n=1)[source]
This performs the same update as
NeighborhoodUpdateSlow
but is streamlined to eliminate calls, to calculate the change in action directly, and to avoid data movement.Proposals are drawn according to
\[\begin{split}\begin{align} \Delta\phi_x &\sim \text{uniform}(-\texttt{interval_phi}, +\texttt{interval_phi}) \\ \Delta n_\ell &\sim W \times [-\texttt{interval_n}, +\texttt{interval_n}] \end{align}\end{split}\]We pick \(\Delta n_\ell\) to be a multiple of the constraint integer \(W\) so that if the adjacent plaquettes satisfy the winding constraint \(dn \equiv 0 \text{ mod }W\) before the update they satisfy it after as well.
See also
On a small 5×5 example this generator yields about three times as many updates per second than
NeighborhoodUpdateSlow
on my machine. This ratio should improve for larger lattices because the change in action is computed directly and is of fixed cost, rather than scaling with the volume.
One reason the NeighborhoodUpdate suffers is that offering updates to 5 fields at onces (\(\phi\) and the 4 ns that touch it) means that the change in action can be large. We can decouple these proposals.
- class supervillain.generator.villain.SiteUpdate(action, interval_phi=3.141592653589793)[source]
This performs the same update to \(\phi\) as
NeighborhoodUpdate
but leaves \(n\) untouched.Proposals are drawn according to
\[\Delta \phi_x \sim \text{uniform}(-\texttt{interval_phi}, +\texttt{interval_phi})\]
- class supervillain.generator.villain.LinkUpdate(action, interval_n=1)[source]
This performs the same update to \(n\) as
NeighborhoodUpdate
but leaves \(\phi\) untouched.Proposals are drawn according to
\[\begin{align} \Delta n_\ell &\sim W \times [-\texttt{interval_n}, +\texttt{interval_n}] \setminus \{0\} \end{align}\]We pick \(\Delta n_\ell\) to be a multiple of the constraint integer \(W\) so that if the adjacent plaquettes satisfy the winding constraint \(dn \equiv 0 \text{ mod }W\) before the update they satisfy it after as well.
Note
You can run
python supervillain/generator/villain/link.py
to compare a pure \(W=1\)NeighborhoodUpdate
ensemble against an ensemble which also doesLinkUpdates
. Note that adding theLinkUpdate
costs essentially 0 time because all the links are done in parallel and there are no python-level for loops.
When \(W=1\) the combination of the SiteUpdate
and LinkUpdate
are ergodic.
But when \(W>1\) the LinkUpdate
only offers changes to \(n\) by multiples of \(W\) to preserve the constraint \(dn = 0 \text{ mod }W\).
For an ergodic algorithm when \(W>1\) we need to offer ways to change \(n\) by 1 (less than \(W\)) while maintaining the constraint.
We need to make closed updates to \(n\), which can be broken up into ExactUpdate
s (which are automatically closed) and HolonomyUpdate
s.
- class supervillain.generator.villain.ExactUpdate(action, interval_z=1)[source]
The
LinkUpdate
only updates n by multiples of \(W\) in order to preserve the constraint \(dn = 0\; (\text{mod } W\)). Another way to preserve the constraint is to update n around a given site in a coordinated way so that \(dn\) is not changed on any of the neighboring plaquettes.One way to accomplish this coordinated change is to start with a zero-form \(z\) and update \(n\) by \(dz\). Then, \(dn = d^2z = 0\).
Proposals are drawn according to
\[\begin{align} z_x &\sim [-\texttt{interval_z}, +\texttt{interval_z}] \setminus \{0\} \end{align}\]
- class supervillain.generator.villain.HolonomyUpdate(action, interval_h=1)[source]
The
ExactUpdate
can change \(n\) by ±1 (even when \(W>1\)), but it does it in a coordinated way—the changes offered are exact, d(a zero form). No combination of exact updates, however, can create a net winding around the torus.This update offers a \(dn\)-preserving update by changing all of the links in a parallel strip around the lattice (see Fig. 2 of Ref. [2]).
Proposals change a whole strip of links simultaneously by
\[\begin{align} h &\sim [-\texttt{interval_h}, +\texttt{interval_h}] \setminus \{0\} \end{align}\]Because the proposal touches a (linearly) extensive number of variables, this update may frequently be rejected.
The combination of the SiteUpdate
, LinkUpdate
, ExactUpdate
, and HolonomyUpdate
is ergodic even when \(W>1\).
But it can be slow to decorrelate.
As mentioned, the HolonomyUpdate
often rejects because it touches a macroscopic number of variables.
A major issue is that the route across the torus is very rigid: it’s just a straight shot.
Smarter worm algorithms can make high-acceptance updates to many variables across the lattice, which can help overcome critical slowing down.
Remember that in the Villain
case we are trying to sample according to
and that we may directly path-integrate out the Lagrange multiplier \(v\) in favor of a constraint
One observable of interest is the Vortex_Vortex
correlation function,
which poses a tricky problem to evaluate, since if we sample configurations of \(Z\) we integrate \(v\) out first and cannot easily compute the observable. Instead, we can think of adding an insertion of the vortex creation and annihilation operators, and we can absorb them into the action before path-integrating out \(v\) and use the insertions to shift the constraint at the insertions
where \(x\) and \(y\) label plaquettes.
The built-in Vortex_Vortex
observable measures this by picking a path between \(x\) and \(y\) on which to change \(n\) to satisfy the new shifted constraint.
Depending on the parameters this inevitably hits an overlap problem, where we must change so many links between \(x\) and \(y\) that the change in action is very big and the correlator is naturally very small, except on rare configurations where it is big.
To address overlap problems of this kind Prokof’ev and Svistunov invented worm algorithms [3] which operate by introducing defects where constraints may be broken, propagating those defects, and when they meet and annihilate the configuration again obeys the constraint. Consider the mixed regular+path integral \(G\) with unspecified normalization \(N\) (that will cancel from all interesting quatities that follow) that integrates over all sectors of possible constraints
A configuration of \(G\) consists of two plaquettes \(h\) and \(t\), \(\phi\in\mathbb{R}\) on sites and \(n\in\mathbb{Z}\) on links which satisfy the \((h, t)\) constraint.
Worm algorithms sample from the larger space of configurations of \(G\) and, when the head \(h\) and tail \(t\) coincide, each \(G\) configuration is a valid \(Z\) configuration appearing with the right relative frequency. This is very different from thinking that we need to carefully design updates to always maintain the constraint. No! Violate the constraint as part of the evolution!
If we want to use a worm as an update in our Markov chain, we need to emit configurations that satisfy the original winding constraint. We can imagine constructing a Markov chain to sample configurations of both \(Z\) and \(G\). To go from a \(Z\) configuration insert the head and tail on the same randomly-chosen location (in this case, a plaquette); the change in action is 0 and this is automatically accepted. This configuration now ‘has a worm’ even though the fields have not changed.
Next we evolve the configuration in the larger set of all \(G\) configurations which need not satisfy the original winding constraint but instead violates it in exactly two places (the head and the tail), in opposite ways. Starting from a ‘diagonal’ configuration which has \(h=t\), we can move the head by one plaquette. As the head moves it changes the \(n\) on the link it crosses so that it changes \(dn\) on both adjacent plaquettes. With a clever choice we can ensure that as the head leaves a plaquette it restores the constraint there and breaks it on the destination plaquette. In this way the head moves around the lattice. This change of \(n\) changes the action \(S\), and so the moves need to be Metropolis-tested in some sense.
Finally, when the head reaches the tail, the constraint is restored everywhere. We might allow the worm to keep evolving, but since it satisfies the winding constraint it is possible to reach a \(Z\) configuration too. To go from a \(G\) configuration to a \(Z\) configuration we require the head and tail to coincide and satisfy the winding constraint; the change in action is 0 and the acceptance is likewise automatically accepted if this change is proposed. Now we have a \(Z\) configuration and add it to our Markov chain (or update it with different generators first). Inside the worm algorithm we need not provide a direct route from any \(Z\) configuration to any other; the only way between \(Z\) configurations is through \(G\), though nothing prevents us from using other generators to go directly from \(Z\) configuration to \(Z\) configuration.
Notice that
where the expectation values are with respect to configurations drawn from \(G\) (not \(Z\)!).
Let us understand this expression for the vortex correlator.
It says that if we draw from the larger space of \(G\) configurations and make a histogram in \(x\) and \(y\), we can normalize that histogram by its value at zero displacement to get the \(V_{x,y}\).
We can measure the histogram as we go and report with the \(\phi\) and \(n\).
The reason to save the histogram as a sample rather than to just accumulate a histogram for the whole worm’s evolution is that once we reach a \(Z\) configuration we update the other variables; that histogram is conditional on those variables; when they change the histogram will too.
So, the worm’s displacement histogram can be saved inline as Vortex_Vortex
, as long as we remember to normalize any DerivedQuantity
that depends on it by the element of the expectation at the origin.
- class supervillain.generator.villain.worm.ClassicWorm(S)[source]
This implements the classic worm of Prokof’ev and Svistunov [3] for the Villain links \(n\in\mathbb{Z}\) which satisfy \(dn \equiv 0 \) (mod W) on every plaquette.
On top of a constraint-satisfying configuration we put down a worm and let the head move, changing the crossed links. We uniformly propose a move in all 4 directions and Metropolize the change.
Additionally, when the head and tail coincide, we allow a fifth possible move, where we remove the worm and emit the updated \(z\) configuration into the Markov chain.
As we evolve the worm we tally the histogram that yields the
Vortex_Vortex
correlation function.Warning
This update algorithm is not ergodic on its own. It doesn’t change \(\phi\) at all and even leaves \(dn\) alone (while changing \(n\) itself). It can be used, for example,
Sequentially
with theSiteUpdate
andLinkUpdate
for an ergodic method.Note
Because it doesn’t change \(dn\) at all, this algorithm can play an important role in sampling the \(W=\infty\) sector, where all vortices are completely killed, though updates to \(\phi\) would still be needed.
Note
This class contains kernels accelerated using numba.
See also
There is
a reference implementation without any numba acceleration
.- inline_observables(steps)[source]
The worm algorithm can measure the
Vortex_Vortex
correlator. We also store theWorm_Length
for each step.
The worm is not ergodic on its own—it doesn’t update \(\phi\), for example, and it cannot change a link by ±W.
But, in combination of with SiteUpdate
and LinkUpdate
it is ergodic;
the worm can replace the combination of ExactUpdate
and HolonomyUpdate
.
The ExactUpdate
can be understood as a very simple worm that takes the tightest nontrivial path,
while the HolonomyUpdate
can be understood as a worm that goes once around the world.
The worm offers dynamically determined constraint-preserving updates and is much more flexible.
In can change the holonomy, for example, by finding a route around the torus that isn’t a straight shot but
runs through the valley of the action.
Finally, we provide a convenience function which provides an ergodic generator.
- supervillain.generator.villain.Hammer(S, worms=1)[source]
The Hammer is just syntactic sugar for a
Sequentially
applied ergodic combination of generators. It may change from version to version as new generators become available or get improved.Note
When \(W=\infty\) we only include updates that leave \(dn=0\) (NOT \(\text{mod }W\)!).
- Parameters
S (a Villain action)
worms (int) – A positive integer saying how many worms to do per iteration.
- Return type
An ergodic generator for updating Villain configurations.
The Worldline Formulation
In the Worldline
formulation the constraint \(\delta m = 0\) restricts which kinds of updates we could propose.
For example, changing only a single link is guaranteed to break the constraint on both ends.
So, we need clever generators to maintain the constraint.
- class supervillain.generator.worldline.PlaquetteUpdate(action)[source]
Ref. [4] suggests a simple update scheme where the links surrounding a single plaquette are updated in concert so that the
Worldline
constraint is maintained.The links on a plaquette are updated in the same (oriented) way, which guarantees that \(\delta m\) is invariant before and after the update, so that if it started 0 everywhere so it remains. The coordinated change is randomly chosen from ±1.
Warning
HOWEVER this algorithm is not ergodic on its own. The issue is that no proposal can change the worldline
TorusWrapping
. Instead, if you start cold with \(m=0\), which has global wrapping of (0,0) you stay in the (0,0) sector.
We can decouple the proposals in the PlaquetteUpdate
, and update the vortex fields and the worldlines independently.
- class supervillain.generator.worldline.VortexUpdate(action, interval_v=1)[source]
This performs the same update to \(v\) as
PlaquetteUpdate
but leaves \(m\) untouched.Proposals are drawn according to
\[\begin{split}\begin{align} \Delta v_p &\sim [-\texttt{interval_v}, +\texttt{interval_v}] \setminus \{0\} &&(W<\infty) \\ \Delta v_p &\sim \text{uniform}(-\texttt{interval_v}, +\texttt{interval_v}) &&(W=\infty) \end{align}\end{split}\]on each plaquette \(p\) independently.
Warning
This update is not ergodic on its own, since it does not change \(m\) at all.
- class supervillain.generator.worldline.CoexactUpdate(action, interval_t=1)[source]
One way to guarantee that \(\delta m = 0\) is to change \(m\) by \(\delta t\) where \(t\) is an integer-valued two-form.
Proposals are drawn according to
\[\begin{align} t_p &\sim [-\texttt{interval_t}, +\texttt{interval_t}] \setminus \{0\} & \Delta m_\ell &= (\delta t)_\ell \end{align}\]Warning
This algorithm is not ergodic on its own. It does not change \(v\) (see the
VortexUpdate
) nor can it produce all coclosed changes—only coexact changes. For a coïnexact coclosed update consider theWrappingUpdate
or theworm
.
To have a fully ergodic algorithm we will also need to update the TorusWrapping
of the worldlines.
- class supervillain.generator.worldline.WrappingUpdate(action, interval_w=1)[source]
Because
CoexactUpdate
fails to change the wrapping, we should separately offer wrapping-changing proposals.We propose coordinated changes on all the x-direction links on a single timeslice and coordinated changes on all the t-direction links on a single spatial slice.
That is, on all the \(m\) on a cycle around the torus we propose to change \(m\) according to
\[\Delta m \sim [-\texttt{interval_w}, +\texttt{interval_w}] \setminus \{0\}\]Warning
This algorithm is not ergodic on its own. The issue is that no proposal can generate coexact changes.
The combination of the VortexUpdate
, CoexactUpdate
, and WrappingUpdate
are ergodic.
However, the may be suboptimal. In particular the WrappingUpdate
touches a large number of links which often leads to very large changes in action
and rejection. Just as in the Villain case we can make smarter updates to a dynamically-determined set of variables with high acceptance by using a worm algorithm.
Unlike the Villain formulation, the Worldline formulation has a constraint even when \(W=1\), \(\delta m = 0\) everywhere, from path-integrating \(\phi\)
Thinking about the Spin_Spin
correlation function, we want to insert \(e^{i(\phi_x-\phi_y)}\), which shifts the constraint to
and we can perform updates in the mixed regular+path integral \(G\) that integrates over all sectors of possible constraints
As in the Villain case, when the head \(h\) and tail \(t\) coincide a configuration of \(G\) is a valid \(Z\) configuration and each appears with the correct relative frequency. Again, the philosophy is to evolve in a larger space with the constraint lifted and celebrate when we receive a constraint-satisfying configuration.
Just like the Villain case we can measure a two-point correlator inline, but in this case the constraint is \(\delta m = 0\) everywhere and constraint-violating insertions are of \(e^{\pm i \phi}\) (as in Spin_Spin
),
which amounts to constructing a normalized histogram after ensemble averaging.
- class supervillain.generator.worldline.worm.ClassicWorm(S)[source]
This implements the classic worm of Prokof’ev and Svistunov [3] for the worldline links \(m\in\mathbb{Z}\) which satisfy \(\delta m = 0\) on every site.
On top of a constraint-satisfying configuration we put down a worm and let the head move, changing the crossed links. We uniformly propose a move in all 4 directions and Metropolize the change.
Additionally, when the head and tail coincide, we allow a fifth possible move, where we remove the worm and emit the updated \(z\) configuration into the Markov chain.
As we evolve the worm we tally the histogram that yields the
Spin_Spin
correlation function.Warning
When \(W>1\) this update algorithm is not ergodic on its own. It doesn’t change \(v\) at all. However, when \(W=1\) we can always pick \(v=0\) (any other choice may be absorbed into \(m\)), and this generator can stand alone.
Note
This class contains kernels accelerated using numba.
See also
There is
a reference implementation without any numba acceleration
.- inline_observables(steps)[source]
The worm algorithm can measure the
Spin_Spin
correlator. We also store theWorm_Length
for each step.
Finally, we provide a convenience function which provides an ergodic generator.
- supervillain.generator.worldline.Hammer(S, worms=1)[source]
The Hammer is just syntactic sugar for a
Sequentially
applied ergodic combination of generators. It may change from version to version as new generators become available or get improved.- Parameters
S (a Worldline action)
worms (int) – A positive integer saying how many worms to do per iteration.
- Return type
An ergodic generator for updating Villain configurations.
Combining Generators
You can combine generators. One simple combination is just the sequential application.
- class supervillain.generator.combining.Sequentially(generators)[source]
Sequentially applies the list of generators as a single step.
For example we can get an ergodic
Worldline
update scheme by combining thePlaquetteUpdate
andWrappingUpdate
.>>> p = supervillain.generator.worldline.PlaquetteUpdate(S) >>> h = supervillain.generator.worldline.WrappingUpdate(S) >>> g = supervillain.generator.combining.Sequentially((p, h))
- Parameters
generators (iterable of generators) – The ordered iterable of generators to apply one after the next.
- step(cfg)[source]
Apply each generator’s
step
one after the next and return the final configuration.
Another (trivial) combination is a decorrelating application.
- class supervillain.generator.combining.KeepEvery(n, generator, blocked_inline=True)[source]
To decorrelate and get honest error estimates it can be helpful to do MCMC but then only analyze evenly-spaced configurations. Rather than keep every single generated configuration and then throw a bunch away, we can keep only those we might analyze.
Note
The number of updates per second will decrease by a factor of n, but the autocorrelation time should be n less. Generating a fixed number of configurations will take n times longer.
>>> p = supervillain.generator.worldline.PlaquetteUpdate(S) >>> g = supervillain.generator.combining.KeepEvery(10, p)
- Parameters
n (int) – How many generator updates to make before emitting a configuration.
generator – The generator to use for updates.
blocked_inline (bool) – Rather than just keeping the
inline_observables()
from the last update, all of the inline measurements are averaged across all n updates. This helps capture rare-but-important measurements that would otherwise be missed.
Monitors
We can make pass-through generators to perform specific tasks, like logging.
- class supervillain.generator.monitor.Logger(generator, channel)[source]
Mostly good for debugging.
import logging logger = logging.getLogger(__name__) ... L = supervillain.lattice.Lattice2D(N) S = supervillain.action.Worldline(L, kappa, W) G = supervillain.generator.worldline.Hammer(S, L.sites) M = supervillain.generator.monitor.Logger(G, logger.info) E = supervillain.Ensemble(S).generate(100, M)
- Parameters
generator (any generator) – The generator to wrap.
channel (some function) – Will be applied to the resulting configuration.
One could construct monitors that, for example, measured observables inline.