Population-dynamics method with a multicanonical feedback control

We discuss the Giardinà-Kurchan-Peliti population dynamics method for evaluating large deviations of time-averaged quantities in Markov processes [Phys. Rev. Lett. 96, 120603 (2006)]. This method exhibits systematic errors which can be large in some circumstances, particularly for systems with weak noise, with many degrees of freedom, or close to dynamical phase transitions. We show how these errors can be mitigated by introducing control forces within the algorithm. These forces are determined by an iteration-and-feedback scheme, inspired by multicanonical methods in equilibrium sampling. We demonstrate substantially improved results in a simple model, and we discuss potential applications to more complex systems.

Figure 1

(a) Trajectories $x^a\left(t\right)$ generated by population dynamics at fixed total population $N_c=4$ for the model system described in Sec. 2c ($\epsilon =1,h=1$). The different colors or line types, the green (dark gray) solid line, green dashed line, yellow (light gray) solid line, and yellow dashed line, represent different copies. At certain times, some copies of the system are removed ($\times$ symbols) and others are duplicated ($\circ$ symbols). The time interval $\mathrm{\Delta }T$ for the cloning procedure is set to be 0.05, and the time step for solving the Langevin equation is 0.001 (see Appendix pp1-s1 for the details of the algorithm). (b) Representative sample paths ${\tilde{x}}^a\left(t\right)$ for the distribution $P_h\left[X\right]$, derived from those in (a) by keeping only trajectories surviving up to final time $\tau =30$. For each cloning event, we also copy the history of the trajectory, which replaces the history of the eliminated trajectory. This means that the trajectories ${\left({\tilde{x}}^a\left(t\right)\right)}_{a=1}^{N_c}$ overlap, especially for early times. For example, in panel (b), the point A appears in the past of the four points B1, $\cdots$, B4. For any point $x^a\left(t\right)$ (such as A, B1, B2, $\cdots$), we define the multiplicity $m_a(t,\tau )$ as the number of trajectories that include this point and survive until the final time $\tau$. So for point $A$, the multiplicity is $m_a(t,\tau )=4$, but for B1, $\cdots$, B4 then $m_a(t,\tau )=1$. [For all points in the trajectories that die before $\tau$, which are not drawn in panel (b), $m_a(t,\tau )=0$.]

Figure 2

(a–d) Distributions $p_{\mathrm{end}}\left(x\right)$ and $p_{\mathrm{ave}}\left(x\right)$, defined in (8) and (9), calculated from the population dynamics method, with various numbers of clones $N_c$. The different panels correspond to a different value of $h$ ($h=\pm 1$) or a different distribution function [$p_{\mathrm{end}}\left(x\right)$ or $p_{\mathrm{ave}}\left(x\right)$]: (a) $p_{\mathrm{end}}\left(x\right)$ for $h=-1$, (b) $p_{\mathrm{ave}}\left(x\right)$ for $h=-1$, (c) $p_{\mathrm{end}}\left(x\right)$ for $h=1$, and (d) $p_{\mathrm{ave}}\left(x\right)$ for $h=1$. For all panels, we set $\epsilon =1$ and $\tau =30$. The numerically exact result is plotted as a black line. We repeat the simulation $1200/N_c$ times, and the result is the average of them (this procedure means that we vary $N_c$ while keeping a fixed computational cost). The results of the population dynamics converge to the analytical ones as $N_c$ increases. (e) $p_{\mathrm{ave}}\left(x\right)$ for $h=1$ (improved estimation) calculated from a population dynamics method with control-with-feedback, as described in Secs. 4c and 5. Results are shown after two iterations of the feedback procedure. The exact distribution $p_{\mathrm{ave}}\left(x\right)$ is again shown as a black line. The comparison between (d) and (e) indicates that the convergence with respect to $N_c$ is improved significantly by the control-with-feedback method. The variance $m_2$ and the relative entropy $D_2$ defined in (12) and (13) both measure how much large values of $N_c$ are required for the cloning procedure to be reliable. For panel (b), (d) and (e), these values are ($m_2=0.068, D_2=0.039$), ($m_2=0.33, D_2=0.17$), and ($m_2=0.0064, D_2=0.0032$) respectively.

Figure 3

The number of independent (distinct) clones ${\tilde{N}}_c\left(t\right)$ obtained from the normal population dynamics method [green line (dark gray) line] for $h=-1$ (a) and $h=1$ (b). The line type corresponds to the value of noise intensity: $\epsilon =1$ (solid line) and $\epsilon =0.1$ (dashed line). We set $N_c=20$ and $\tau =30$. When the distributions $p_{\mathrm{ave}}$ and $p_{\mathrm{end}}$ are very different from each other, we expect that ${\tilde{N}}_c\left(t\right)$ decreases rapidly as $t$ decreases from $\tau$: to illustrate this, note that (for $\epsilon =1$) $m_2=0.068$ for $h=-1$ and $m_2=0.33$ for $h=1$: the same ordering is preserved for smaller $\epsilon$. We also plot ${\tilde{N}}_c\left(t\right)$ obtained from the controlled population dynamics [yellow line (light gray) line] with the control-with-feedback explained in Secs. 4c and 5. The larger values of ${\tilde{N}}_c\left(t\right)$ obtained with the control-with-feedback method lead to smaller statistical uncertainties in the results.

Figure 4

Estimates of $\tilde{G}(\tilde{h}=1)$, as $\epsilon$ is varied. We compare results from the normal population dynamics and from the control-with-feedback method explained in Secs. 4c and 5. The analytical result for ${lim}_{\epsilon \to 0}\tilde{G}\left(\tilde{h}\right)$ is shown as a red dashed line, and the characteristic value of the noise intensity ${\epsilon }^{*}$, defined in (14), is plotted as a green vertical solid line. The standard method fails for $\epsilon$ smaller than ${\epsilon }^{*}$, but the control-with-feedback method (black continuous line and black circles) converges to the correct value even for $\epsilon <{\epsilon }^{*}$.

Figure 5

A schematic map illustrating the methodological situation of the controlled population dynamics.

Figure 6

The integrands of $m_2$ and $D_2$ defined as (a) $p_{\mathrm{end}}^w\left(x\right)\left\{{\left[\frac{p_{\mathrm{ave}}^w\left(x\right)}{p_{\mathrm{end}}^w\left(x\right)}\right]}^2-1\right\}$ and (b) $p_{\mathrm{ave}}^w\left(x\right)log\left[\frac{p_{\mathrm{ave}}^w\left(x\right)}{p_{\mathrm{end}}^w\left(x\right)}\right]$ [see (12) and (13)] for the standard (normal) population dynamics method ($w=0$) and for the control-with-feedback method ($w$: obtained from the control-with-feedback method). For the control-with-feedback method, we set $N_c=20$. In the legends, we put the values corresponding to $m_2$ and $D_2$.

Figure 7

Plots of the polynomial $3x^5-4hx-2h$ for several $h$. The roots of this polynomial determine the concentration points of $p_{\mathrm{ave}}\left(x\right)$ for $\epsilon \to 0$ in the model system considered.

Figure 8

(a) Generating function ${lim}_{\epsilon \to 0}{G}_{\epsilon }\left(h\right)$. (b) The first derivatives, where the green (dark gray) solid line represents $(d/dh){lim}_{\epsilon \to 0}{G}_{\epsilon }\left(h\right)$, the red dotted line represents $-(d/dh){lim}_{\epsilon \to 0}{G}_{\epsilon }(-h)$, and the black solid line represents a straight line $h^{1/5}$. We find that the first derivative converges to 0 as a power law $\sim h^{1/5}$ (as can also be checked analytically). (c) The second derivative of ${lim}_{\epsilon \to 0}{G}_{\epsilon }\left(h\right)$ with respect to $h$. These are calculated from (D11, D12). We can see that the second derivative shows a singularity at $h=0$, although the first derivative converges to 0. This represents a second-order dynamical phase transition.

Figure 9

The functions $\partial W_{\mathrm{F}}(x,t)/\partial x$ (a, b) and $\partial W_{\mathrm{B}}(x,t)/\partial x$ (c, d) obtained in the large $t$ limit by solving numerically (D15) and (D16) [yellow lines (light gray) lines]. We set $h=1$ (a, c) and $h=-1$ (b, d). The line types correspond different values of $\epsilon$: dash-dotted, dashed, and solid lines correspond to $\epsilon =1,0.5,0.1$, respectively. To illustrate the determination of the $\pm$ sign of $C_h$ in the analytical results (D17) and (D18), we also plot on each subfigure those results with the choice of $C_h=1$ (for all $x$) as black solid lines and the choice of $C_h=-1$ (for all $x$) as black dashed lines. As the noise goes to zero, we observe the convergence of the functions $\partial W_{\mathrm{F}}(x,t)/\partial x$ and $\partial W_{\mathrm{B}}(x,t)/\partial x$ determined numerically at large $t$ towards the analytical line (D17) and (D18), where the $+$ sign in $\pm$ is taken for $x<x_{\min }$ and the $-$ sign is taken for $x>x_{\min }$.