# Dynamical Lee–Yang zeros for continuous-time and discrete-time stochastic
processes

Hiroki Yoshida  Department of Physics, Tokyo Institute of Technology, Tokyo
152–8551, Japan  Kazutaka Takahashi  Institute of Innovative Research, Tokyo
Institute of Technology, Kanagawa 226–8503, Japan

###### Abstract

We describe classical stochastic processes by using dynamical Lee–Yang zeros.
The system is in contact with external leads and the time evolution is
described by the two-state classical master equation. The cumulant generating
function is written in a factorized form and the current distribution is
characterized by the dynamical Lee–Yang zeros. We show that a continuous
distribution of zeros is obtained by discretizing the time variable. When the
transition probability is a periodically oscillating function of time, the
distribution of zeros splits into many parts. We study the geometric property
of the current by comparing the result with that of the adiabatic
approximation. We also use the Floquet–Magnus expansion in the continuous-time
case to study dynamical effects on the current at the fast-driving regime.

##  I Introduction

Understanding statistical fluctuations of current is one of the main
objectives in nonequilibrium physics. When the system is coupled to multiple
reservoirs, we observe a current through the system. The fluctuation theorem
tells us how the underlying symmetry of the system is reflected to the current
distribution [1, 2, 3].

In the long-time limit, the system settles down to a stationary behavior
irrespective of the choice of the initial condition. Full counting statistics
describes distributions of transferred charges through the system [2, 4, 5, 6,
7]. The generating function can be treated as an analog of the partition
function in equilibrium statistical mechanics.

The analytic properties of the partition function can be described by the
Lee–Yang zeros [8, 9, 10]. The partition function is a positive quantity and
goes to zero only when we set the parameters to unphysical values. The only
exception is when the system has a phase transition. It involves a breakdown
of the analyticity and can be described by the distribution of the zeros.

A similar behavior is also expected to hold in the generating function of the
full counting statistics. Previous studies found that the generating function
is described by the dynamical Lee–Yang zeros [11, 12, 13, 14, 15, 16, 17, 18].
We note that the Lee–Yang zeros are useful not only in the presence of the
phase transitions but also for general systems. By knowing the distribution of
zeros, we can obtain the complete information of the system. Any statistical
quantities are represented by zeros. It was shown in quantum systems [19, 20]
and in classical stochastic systems [17] that the Lee–Yang zeros are not
artificial mathematical objects but are directly related to the experimental
observables.

Zeros of the generating function can be clearly understood when the number of
events is finite, or countably infinite. For a stochastic process described by
the probability distribution $\\{P_{n}\\}_{n=-M,-M+1,\dots,M}$, the generating
function is defined as

$\displaystyle Z(\chi)=\sum_{n=-M}^{M}(e^{\chi})^{n}P_{n}.$ (1)

When the positive integer parameter $M$ is finite, the generating function,
written as

$\displaystyle
Z(\chi)=\frac{1}{(e^{\chi})^{M}}\prod_{i=1}^{2M}\frac{e^{\chi}-z_{i}}{1-z_{i}},$
(2)

is completely characterized by a set of zeros $\\{z_{i}\\}_{i=1,\dots,2M}$.
Each $z_{i}$ represents a point where the generating function goes to zero in
complex plane of $z=e^{\chi}$. Statistical quantities calculated from the
generating function are represented by using the zeros.

The number of zeros becomes large as the number of events increases. However,
previous studies found that the number of zeros in stochastic systems
described by the two-state classical master equation is given by a small
finite value [12, 13, 15, 17]. This property is totally different from the
original Lee–Yang theory where many zeros are accumulated to form a continuous
distribution in a complex plane in the thermodynamic limit. To find a similar
behavior in stochastic systems, we need to find a controllable parameter that
determines the number of zeros.

A different interesting problem is the situation when the system is driven
periodically. Then, the transition-rate matrix in the master equation is time
dependent. When we treat the system by the adiabatic approximation, the
generating function is defined at each time and the zeros oscillate as a
function of time. However, they do not represent the zeros of the global
generating function. It would be an interesting problem to clarify how to
describe the periodically driven systems by the dynamical Lee–Yang zeros.

In this paper, we study distributions of the dynamical Lee–Yang zeros by using
the two-state classical master equation in several situations. We show that
the number of zeros is controllable by discretizing the time variable. Then,
we can find a nontrivial distribution of zeros in the long-time limit. To
study periodically driven systems, we use the adiabatic approximation in the
slow-driving regime [21, 22] and the Floquet–Magnus expansion in the fast-
driving regime [23, 24].

The organization of the paper is as follows. In Sec. II, we treat the
continuous-time master equation to reproduce known results on the dynamical
Lee–Yang zeros. The analysis is extended to the discrete-time master equation
to obtain a nontrivial distribution of zeros in Sec. III. Then, we study a
periodically driven system from an adiabatic picture in Sec. IV. We also study
the continuous-time system with a high frequency in Sec. V. Section VI is
devoted to summary.

##  II Continuous-time process

Figure 1: Two-state classical stochastic processes. The system takes two
states $|0\rangle$ and $|1\rangle$. Each of $a_{\rm L}$, $b_{\rm L}$, $a_{\rm
R}$, and $b_{\rm R}$ represents a transition rate.

When we pursue an analogy between the current generating function for
stochastic processes and the partition function in statistical mechanics, it
is reasonable to study the long-time behavior of the stochastic processes. The
time evolution of the probability distribution is generally given by the
master equation under the assumption that the process is Markovian. To make
the discussion concrete, we treat the system depicted in Fig. 1. The system
takes two states $|0\rangle=(1,0)^{\rm T}$ (“empty”) and $|1\rangle=(0,1)^{\rm
T}$ (“filled”), and we observe flows among left and right leads. When the
transition-rate matrix in the master equation is time independent, the
dynamical Lee-Yang zeros were obtained in previous studies [12, 13, 15, 17].
In this section, we summarize the basic definitions and notations, and
reproduce the known results together with several small findings.

The time evolution of the probability distribution of states
$|p(t)\rangle=p_{0}(t)|0\rangle+p_{1}(t)|1\rangle$ is described by the master
equation. We introduce the counting field $\chi$ to measure flows between the
system and the right lead and define the modified master equation [22]

$\displaystyle\frac{\partial}{\partial
t}|p^{\chi}(t)\rangle=w^{\chi}|p^{\chi}(t)\rangle,$ (3)

with the transition-rate matrix

$\displaystyle w^{\chi}=\left(\begin{array}[]{cc}-b&a^{\chi}\\\
b^{\chi}&-a\end{array}\right)=\left(\begin{array}[]{cc}-(b_{\rm L}+b_{\rm
R})&a_{\rm L}+a_{\rm R}e^{\chi}\\\ b_{\rm L}+b_{\rm R}e^{-\chi}&-(a_{\rm
L}+a_{\rm R})\end{array}\right),$ (8)

where each of $a_{\rm L}$, $b_{\rm L}$, $a_{\rm R}$, and $b_{\rm R}$
represents a transition rate and takes a positive value. The generating
function is obtained from $|p^{\chi}(t)\rangle$ as $\langle
I|p^{\chi}(t)\rangle$ where $\langle I|=\sum_{n=0}^{1}\langle n|$. Since the
probabilistic properties of $|p(t)\rangle$ are maintained throughout the time
evolution, the off-diagonal components of $w^{\chi=0}$ are non-negative and
the diagonal components are determined from the relation $\langle
I|w^{\chi=0}=0$. The counting field is incorporated only in the off-diagonal
components. To describe the long-time limit, we define the cumulant generating
function

$\displaystyle g(\chi)=\lim_{t\to\infty}\frac{1}{t}\ln\langle
I|p^{\chi}(t)\rangle.$ (9)

The $k$th moment of the current is given by

$\displaystyle
J_{k}=\left.\frac{\partial^{k}}{\partial\chi^{k}}g(\chi)\right|_{\chi=0}.$
(10)

The current generating function is obtained from the Legendre transformation
of $g(\chi)$ as we demonstrate below.

The modified master equation is solved by using the spectral representation of
$w^{\chi}$ as

$\displaystyle
w^{\chi}=\sum_{n=0}^{1}\lambda_{n}^{\chi}|R_{n}^{\chi}\rangle\langle
L_{n}^{\chi}|,$ (11)

where $\lambda_{n}^{\chi}$ denotes the eigenvalue, and $\langle L_{n}^{\chi}|$
and $|R_{n}^{\chi}\rangle$ are the corresponding left and right eigenstates
respectively. The eigenstates satisfy the orthonormal relation $\langle
L_{m}^{\chi}|R_{n}^{\chi}\rangle=\delta_{m,n}$ and the resolution of unity
$\sum_{n}|R_{n}^{\chi}\rangle\langle L_{n}^{\chi}|=1$. Although the
diagonalization of the general transition-rate matrix is not always possible,
we can find the spectral representation in the present two-state case as we
see below. The eigenvalues are explicitly written as

$\displaystyle\lambda_{0,1}^{\chi}=-\frac{a+b}{2}\pm\frac{a+b}{2}\sqrt{\frac{(e^{\chi}-z_{1})(e^{\chi}-z_{2})}{e^{\chi}(1-z_{1})(1-z_{2})}},$
(12)

where $z_{1}$ and $z_{2}$ are given by

$\displaystyle z_{1,2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left[\frac{(a+b)^{2}}{4b_{\rm L}a_{\rm
R}}-1-\frac{a_{\rm L}b_{\rm R}}{b_{\rm L}a_{\rm R}}\right]$ (13)
$\displaystyle\pm\sqrt{\frac{1}{4}\left[\frac{(a+b)^{2}}{4b_{\rm L}a_{\rm
R}}-1-\frac{a_{\rm L}b_{\rm R}}{b_{\rm L}a_{\rm R}}\right]^{2}-\frac{a_{\rm
L}b_{\rm R}}{b_{\rm L}a_{\rm R}}}.$

The modified master equation for a given initial condition $|p(0)\rangle$ is
solved to give the generating function as

$\displaystyle\langle I|p^{\chi}(t)\rangle=e^{\lambda_{0}^{\chi}t}\langle
I|R_{0}^{\chi}\rangle\langle
L_{0}^{\chi}|p(0)\rangle+e^{\lambda_{1}^{\chi}t}\langle
I|R_{1}^{\chi}\rangle\langle L_{1}^{\chi}|p(0)\rangle.$

The component $n=0$ represents the stationary state since
$\lambda_{0}^{\chi}>\lambda_{1}^{\chi}$ holds for arbitrary real values of
$\chi$. We also note that $\lambda_{0}^{\chi=0}=0$ and
$\lambda_{1}^{\chi=0}=-(a+b)<0$.

We are interested in zeros of the generating function by extending $\chi$ to
complex values. They are obtained by solving

$\displaystyle\lambda_{0}^{\chi}-\lambda_{1}^{\chi}=\frac{1}{t}\ln\left(-\frac{\langle
I|R_{1}^{\chi}\rangle\langle L_{1}^{\chi}|p(0)\rangle}{\langle
I|R_{0}^{\chi}\rangle\langle L_{0}^{\chi}|p(0)\rangle}\right).$ (15)

Although the general solution is dependent on $t$ and on the initial
condition, the equation reduces to $\lambda_{0}^{\chi}=\lambda_{1}^{\chi}$ at
large $t$. From the explicit forms of the eigenvalues in Eq. (12), we conclude
that the dynamical Lee–Yang zeros are given by $e^{\chi}=z_{1}$ and
$e^{\chi}=z_{2}$. These points are negative, which means that the zeros are
located on inaccessible region in complex plane of $z=e^{\chi}$. The cumulant
generating function is given by

$\displaystyle
g(\chi)=\lambda_{0}^{\chi}=-\frac{a+b}{2}+\frac{a+b}{2}\sqrt{\frac{(e^{\chi}-z_{1})(e^{\chi}-z_{2})}{e^{\chi}(1-z_{1})(1-z_{2})}}.$
(16)

This expression was obtained in previous works [12, 13, 15, 17].

The representation of the cumulant generating function by zeros is convenient
to find nontrivial relations of the current distributions. We write $J_{k}$ in
Eq. (10) as $J_{k}=(a+b)j_{k}/2$. Then, $j_{k}$ is only dependent on $z_{1}$
and $z_{2}$. For example, we have

$\displaystyle
j_{1}=\frac{1}{2}\left(\frac{1}{1-z_{1}}+\frac{1}{1-z_{2}}-1\right).$ (17)

The higher-order moments are calculated in a similar way. As they are
parametrized only by $z_{1}$ and $z_{2}$, we can find nontrivial relations
between the moments. Up to the fourth order of the cumulants, we obtain

$\displaystyle j_{3}+3j_{2}j_{1}-j_{1}=0,$ (18) $\displaystyle
j_{4}+3j_{2}^{2}-12j_{2}j_{1}^{2}-j_{2}+3j_{1}^{2}=0.$ (19)

As far as we understand, these relations have not been obtained in previous
studies.

The current distribution function $P(J,t)$ is calculated from the Legendre
transformation $\varphi(J)=g(\chi)-J\chi$ as $P(J,t)\sim\exp(t\varphi(J))$ at
large $t$. The cumulant generating function in Eq. (16) is invariant under the
transformation $e^{\chi}\to z_{1}z_{2}e^{-\chi}$:

$\displaystyle g(-\chi+\ln(z_{1}z_{2}))=g(\chi).$ (20)

In the Legendre transformation, the relation between $J$ and $\chi$ is given
by $J(\chi)=\partial_{\chi}g(\chi)$. Then, $J(-\chi+\ln(z_{1}z_{2}))=-J(\chi)$
and we have

$\displaystyle\varphi(J)-\varphi(-J)=J\ln\sqrt{z_{1}z_{2}}.$ (21)

In the context of nonequilibrium thermodynamics, this relation is known as the
fluctuation theorem. The product of the zeros $z_{1}z_{2}=a_{\rm L}b_{\rm
R}/b_{\rm L}a_{\rm R}$ gives the affinity ${\cal A}=-\ln(z_{1}z_{2})$ which
represents a bias between the left and right leads.

##  III Discrete-time process

The polynomial representation of the generating function in Eq. (2) is
possible only when the number of events is finite. Equation (16) does not have
a simple polynomial form and it is not clear why the number of zeros is given
by 2. We note that the number of zeros is not related to the number of states
in the master equation, as we confirm below.

Here we treat a discrete-time process. It is suitable to describe rare events
where the charge transfer occurs sporadically. We set a finite-time interval
$\Delta t$ and the total process time is discretized as $t=2M\Delta t$. In
this setting, the charge transfer is counted discretely and the number of
events is controllable by the integer parameter $M$. In the following
analysis, we take the long-time limit $M\to\infty$ while keeping $\Delta t$
finite.

The discrete master equation is given by

$\displaystyle|p(k)\rangle=(1+W)|p(k-1)\rangle,$ (22)

where the integer index $k$ denotes a discretized time $t/\Delta t$. Each
element of $1+W$ represents a probability and must be smaller than unity. The
transition-rate matrix $w$ in the continuous limit is obtained by using the
relation $W=w\Delta t$. When we include the counting field $\chi$ to $W$,
$W^{\chi}$ is a linear combination of $e^{\chi}$ and $e^{-\chi}$. We write the
transition matrix

$\displaystyle W^{\chi}=\left(\begin{array}[]{cc}-B&A^{\chi}\\\
B^{\chi}&-A\end{array}\right)=\left(\begin{array}[]{cc}-(B_{\rm L}+B_{\rm
R})&A_{\rm L}+A_{\rm R}e^{\chi}\\\ B_{\rm L}+B_{\rm R}e^{-\chi}&-(A_{\rm
L}+A_{\rm R})\end{array}\right).$ (27) (28)

Then, we see that the generating function $\langle
I|p^{\chi}(2M)\rangle=\langle I|(1+W^{\chi})^{2M}|p(0)\rangle$ can be
factorized as

$\displaystyle\langle
I|p^{\chi}(2M)\rangle=\prod_{k=0}^{M-1}\frac{(e^{\chi}-z_{1}(k))(e^{\chi}-z_{2}(k))}{e^{\chi}(1-z_{1}(k))(1-z_{2}(k))}.$
(29)

As in the continuous-time case, we write $W^{\chi}$ as

$\displaystyle
W^{\chi}=\sum_{n=0}^{1}\Lambda_{n}^{\chi}|R_{n}^{\chi}\rangle\langle
L_{n}^{\chi}|,$ (30)

and the generating function

$\displaystyle\langle I|p^{\chi}(2M)\rangle$ $\displaystyle=$
$\displaystyle(1+\Lambda_{0}^{\chi})^{2M}\langle I|R_{0}^{\chi}\rangle\langle
L_{0}^{\chi}|p(0)\rangle$ (31)
$\displaystyle+(1+\Lambda_{1}^{\chi})^{2M}\langle I|R_{1}^{\chi}\rangle\langle
L_{1}^{\chi}|p(0)\rangle.$

We want to find the zeros of the generating function, $z_{1}(k)$ and
$z_{2}(k)$ in Eq. (29). They are obtained by solving $\langle
I|p^{\chi}(2M)\rangle=0$. In the continuous-time case, the solution is
expressed as in Eq. (15) and the right hand side of the equation is neglected
in the large-time limit. A similar analysis is applied for the present
discrete-time case. A notable difference in that case is that we need to
introduce the factor $(e^{i\pi k/M})^{2M}=1$ with $k=0,1,\dots,M-1$. We obtain
the relation

$\displaystyle 1+\Lambda_{0}^{\chi}=(1+\Lambda_{1}^{\chi})e^{i\pi k/M}.$ (32)

For a given $k$, there exist two zeros, $z_{1}(k)$ and $z_{2}(k)$. They are
obtained by solving the equation

$\displaystyle\sqrt{\frac{(z-z_{1})(z-z_{2})}{z(1-z_{1})(1-z_{2})}}=\frac{i}{R}\tan\frac{\pi
k}{2M},$ (33)

with respect to $z$. Here, $z_{1,2}$ are given by Eq. (13) with the
replacement $(a_{\rm L},b_{\rm L},a_{\rm R},b_{\rm R})\to(A_{\rm L},B_{\rm
L},A_{\rm R},B_{\rm R})$, and $R=(A+B)/[2-(A+B)]$. The explicit forms of
$z_{1,2}(k)$ are given by

$\displaystyle
z_{1,2}(k)=\frac{z_{1}+z_{2}}{2}-\frac{(1-z_{1})(1-z_{2})\tan^{2}\frac{\pi
k}{2M}}{2R^{2}}$
$\displaystyle\pm\sqrt{\left[\frac{z_{1}+z_{2}}{2}-\frac{(1-z_{1})(1-z_{2})\tan^{2}\frac{\pi
k}{2M}}{2R^{2}}\right]^{2}-z_{1}z_{2}}.$

The zeros are located on the negative real axis of $z=e^{\chi}$ with
$z_{1}\leq z_{1}(k)<0$ and $-\infty<z_{2}(k)\leq z_{2}$. We note that
$z_{1}(0)=z_{1}$ and $z_{2}(0)=z_{2}$ represents the zeros in the continuous-
time case. We also note that the relation

$\displaystyle z_{1}(k)z_{2}(k)=z_{1}z_{2}=\frac{A_{\rm L}B_{\rm R}}{B_{\rm
L}A_{\rm R}}$ (35)

holds for any $k$. This property is due to the symmetry of the transition-rate
matrix, as discussed in the continuous-time case.

At the limit $\Delta t\to 0$, we obtain $R\to 0$ and
$(z_{1}(k),z_{2}(k))_{k\neq 0}\to(0,-\infty)$. Contributions from
$(z_{1}(k),z_{2}(k))$ with $k\neq 0$ do not affect the current since they
satisfy

$\displaystyle\frac{1}{1-z_{1}(k)}+\frac{1}{1-z_{2}(k)}-1=\frac{1-z_{1}z_{2}}{(1-z_{1}(k))(1-z_{2}(k))}\to
0$

at the limit. Then, the system is characterized by two zeros $z_{1}$ and
$z_{2}$. This behavior is consistent with the result in the previous section.
To find the nontrivial form of the cumulant generating function in Eq. (16),
we must carefully take the continuous-time limit $\Delta t\to 0$ and the long-
time limit $M\to\infty$ simultaneously.

We can find a continuous distribution of zeros by taking the limit
$M\to\infty$ while keeping $\Delta t$ finite. At the limit, the density of
zeros defined by

$\displaystyle\rho(z)=\frac{1}{2M}\sum_{k=0}^{M-1}\left(\delta(z-z_{1}(k))+\delta(z-z_{2}(k))\right)$
(37)

is calculated from the number of $k$ within the interval $dz$. We obtain

$\displaystyle\rho(z)$ $\displaystyle=$
$\displaystyle\left|\frac{-i}{2\pi}\frac{\partial}{\partial
z}\ln\frac{1+\Lambda_{0}^{\chi=\ln z}}{1+\Lambda_{1}^{\chi=\ln z}}\right|$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\frac{R\sqrt{\frac{-(z-z_{1})(z-z_{2})}{z(1-z_{1})(1-z_{2})}}}{1+R^{2}\frac{-(z-z_{1})(z-z_{2})}{z(1-z_{1})(1-z_{2})}}\left|\frac{1}{z-z_{1}}+\frac{1}{z-z_{2}}-\frac{1}{z}\right|,$

for $z_{1}<z<0$ and $z<z_{2}$, and $\rho(z)=0$ otherwise. This function is
plotted in Fig. 2. The distribution of zeros is characterized by the edge
points $z_{1}$ and $z_{2}$. These points represent zeros in the continuous-
time limit and the affinity is given by their product as ${\cal
A}=-\ln(z_{1}z_{2})$.

The cumulant generating function is obtained from the density of zeros as

$\displaystyle g(\chi)$ $\displaystyle=$
$\displaystyle\frac{1}{2M}\sum_{k=0}^{M-1}\ln\frac{(e^{\chi}-z_{1}(k))(e^{\chi}-z_{2}(k))}{e^{\chi}(1-z_{1}(k))(1-z_{2}(k))}$
(39) $\displaystyle=$ $\displaystyle\int
dz\,\rho(z)\ln\frac{e^{\chi}-z}{1-z}-\frac{1}{2}\chi.$

The current distribution function $P(J,2M)\sim\exp(2M\varphi(J))$ is also
written by the density of zeros. By using the Legendre transformation, we can
write

$\displaystyle\varphi(J)=\int
dz\,\rho(z)\ln\frac{e^{\chi}-z}{1-z}-\left(J+\frac{1}{2}\right)\chi,$ (40)

where the relation between $J$ and $\chi$ is given by

$\displaystyle J=\int dz\,\rho(z)\frac{e^{\chi}}{e^{\chi}-z}-\frac{1}{2}.$
(41)

$\varphi(J)$ is plotted in Fig. 3.

Figure 2: The density of zeros in Eq. (LABEL:doz) for discrete-time stochastic
processes. (a) We choose the parameters in the transition matrix in Eq. (28)
as $(A_{\rm L},B_{\rm L},A_{\rm R},B_{\rm R})=(0.3,0.2,0.3,0.2)$. The affinity
is given by ${\cal A}=-\ln(z_{1}z_{2})=0$. (b) $(A_{\rm L},B_{\rm L},A_{\rm
R},B_{\rm R})=(0.3,0.2,0.4,0.1)$ and ${\cal A}=-\ln(z_{1}z_{2})=-\ln(3/8)>0$.
Figure 3:  The solid lines denote the rate function $\varphi(J)$ in Eq. (40)
with the density of zeros in Fig. 2. We use the parametrization of $W$ in Fig.
2(a) for the line with “${\cal A}=0$”, and the parametrization in Fig. 2(b)
for “${\cal A}>0$”. The marked points denote results from simulations. We take
$M=2000$ and the number of samples is $10^{5}$.

To confirm that the description by the zeros gives a reasonable result, we
calculate the current distribution from numerical simulations. The state of
the system, $|0\rangle$ or $|1\rangle$, at each time $k$ is denoted by
$n(k)=0$ or 1. We generate a uniform random number $r(k)$ with $0\leq r(k)\leq
1$ at each $k$ to write the state $n(k)$ as

$\displaystyle n(k)$ $\displaystyle=$ $\displaystyle(1-n(k-1))\frac{1+{\rm
sgn}(B-r(k))}{2}$ (42) $\displaystyle+n(k-1)\frac{1-{\rm sgn}(A-r(k))}{2}.$

The net flow to the right lead is calculated as

$\displaystyle j(k)$ $\displaystyle=$ $\displaystyle
j(k-1)-(1-n(k-1))\frac{1+{\rm sgn}(B_{\rm R}-r(k))}{2}$ (43)
$\displaystyle+n(k-1)\frac{1+{\rm sgn}(A_{\rm R}-r(k))}{2},$

with the initial condition $j(0)=0$. Then, the current is given by
$J=j(2M)/2M$. We produce discrete-time sequences of $M=2000$ steps from a
given set of transition probabilities. The current distribution
$P(J,2M)\sim\exp(2M\varphi(J))$ is obtained from $10^{5}$ samples. The result
is plotted in Fig. 3. We see that the result is consistent with that from the
zeros.

##  IV Periodically driven process

When the transition probability fluctuates periodically, we observe a
nontrivial distribution of the current. It is an interesting problem to study
the properties of the system by using the dynamical Lee–Yang zeros. Although
the complete study of the dynamical effects is a difficult problem in general,
we can use the adiabatic approximation when the transition probability slowly
changes as a function of time. It is well known in the adiabatic regime that
the nontrivial geometric effect is observed in the current distributions [21,
22]. In this section, we study the dynamical Lee–Yang zeros in periodically
driven processes. As we mentioned in Sec. I, the dynamical Lee–Yang zeros
cannot be found in the continuous-time case when we use the adiabatic
approximation. This problem does not arise in the discrete-time case. By
extending the analysis of the previous section, we study periodically driven
systems in the discrete-time case.

###  IV.1 Discrete-time formulation

We denote the period of the oscillation by an integer $N$ and write the
probability distribution as

$\displaystyle|p^{\chi}(2MN)\rangle=(U_{N}^{\chi})^{2M}|p(0)\rangle,$ (44)

where

$\displaystyle
U_{N}^{\chi}=(1+W^{\chi}(N))(1+W^{\chi}(N-1))\cdots(1+W^{\chi}(1)).$

The transition matrix at discrete time $k=t/\Delta t$ is denoted by $W(k)$ and
satisfies $W(N+1)=W(1)$. The total time is given by $t=2MN\Delta t$. The
generating function is written as in Eq. (29) with the replacement $M\to MN$.
The number of zeros is given by $2MN$. When we use the spectral decomposition

$\displaystyle
U^{\chi}_{N}=\sum_{n=0}^{1}(1+\Lambda_{n}^{\chi})|R_{n}^{\chi}\rangle\langle
L_{n}^{\chi}|,$ (46)

the zeros of the generating function are given by solving Eq. (32).

We consider the case where the transition matrix with $\chi$ is written in a
form of Eq. (28) at each time. Then, the eigenvalues of $U_{N}^{\chi}$ are
written as

$\displaystyle 1+\Lambda_{n}^{\chi}$ $\displaystyle=$
$\displaystyle\frac{1}{2}{\rm Tr}\,U_{N}^{\chi}\pm\left(1-\frac{1}{2}{\rm
Tr}\,U_{N}^{\chi=0}\right)$ (47)
$\displaystyle\times\sqrt{\frac{(e^{\chi}-z_{1})\cdots(e^{\chi}-z_{2N})}{(e^{\chi})^{N}(1-z_{1})\cdots(1-z_{2N})}},$

to define $\\{z_{1},z_{2},\dots,z_{2N}\\}$. By using this form in Eq. (32), we
can find $2N$ zeros $\\{z_{1}(k),z_{2}(k),\dots,z_{2N}(k)\\}$ at each $k$. The
zeros are located in several separated domains where the function in the
square root in Eq. (47) is negative. The number of the domains is given by
$2\times\lceil{N/2}\rceil$. Each point $z_{i}(0)=z_{i}$ ($i=1,2,\dots,2N$)
denotes an edge of a domain.

The density of zeros in the present case is defined as

$\displaystyle\rho(z)=\frac{1}{2MN}\sum_{k=0}^{M-1}\sum_{i=1}^{2N}\delta(z-z_{i}(k)).$
(48)

When we assume that all of $\\{z_{1},z_{2},\dots,z_{2N}\\}$ take real
(negative) values, by taking the limit $M\to\infty$, we can write

$\displaystyle\rho(z)=\left|\frac{-i}{2\pi N}\frac{\partial}{\partial
z}\ln\frac{1+\Lambda_{0}^{\chi=\ln z}}{1+\Lambda_{1}^{\chi=\ln z}}\right|$
(49)

in the domains of definition. The cumulant generating function is calculated
from the second line of Eq. (39).

###  IV.2 Dynamical current and geometrical current

We examine the relation between the current distribution and the density of
zeros. Before studying the discrete-time system, we summarize the result in
the continuous-time system [22, 25, 26]. When the system is driven
periodically, each of $a_{\rm L}$, $b_{\rm L}$, $a_{\rm R}$, and $b_{\rm R}$
is represented as a function of $\theta=\omega t$ with the period $2\pi$. The
current consists of the dynamical part and the geometrical part: $J=J_{\rm
d}+J_{\rm g}$. The dynamical current is written by using the largest
eigenvalue of the transition-rate matrix as

$\displaystyle J_{\rm d}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\chi}\left.\int_{0}^{2\pi}\frac{d\theta}{2\pi}\,\lambda_{0}^{\chi}(\theta)\right|_{\chi=0}$
(50) $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\frac{d\theta}{2\pi}\,\frac{b_{\rm
L}(\theta)a_{\rm R}(\theta)-a_{\rm L}(\theta)b_{\rm
R}(\theta)}{a(\theta)+b(\theta)}.$

This is independent of the frequency $\omega$. The geometrical part is written
by using the expansion with respect to $\omega$. At the first order, the
geometrical current is written by using the Berry curvature as

$\displaystyle J_{\rm g}$ $\displaystyle=$
$\displaystyle-\omega\frac{\partial}{\partial\chi}\left.\int_{0}^{2\pi}\frac{d\theta}{2\pi}\,\langle
L_{0}^{\chi}(\theta)|\partial_{\theta}R_{0}^{\chi}(\theta)\rangle\right|_{\chi=0}$
(51) $\displaystyle=$
$\displaystyle\omega\int_{0}^{2\pi}\frac{d\theta}{2\pi}\,\frac{a_{\rm
R}(\theta)+b_{\rm
R}(\theta)}{a(\theta)+b(\theta)}\frac{d}{d\theta}\frac{a(\theta)}{a(\theta)+b(\theta)}.$

The current is represented by a flux penetrating a surface in parameter space
[22, 27].

To examine the geometric property of the system, it is useful to write each
parameter as follows:

$\displaystyle\frac{a_{\rm
L}(\theta)}{a+b}=\frac{1+r(\theta)}{2}\sin^{2}\left(\frac{\phi_{a}(\theta)}{2}\right),$
(52) $\displaystyle\frac{a_{\rm
R}(\theta)}{a+b}=\frac{1+r(\theta)}{2}\cos^{2}\left(\frac{\phi_{a}(\theta)}{2}\right),$
(53) $\displaystyle\frac{b_{\rm
L}(\theta)}{a+b}=\frac{1-r(\theta)}{2}\sin^{2}\left(\frac{\phi_{b}(\theta)}{2}\right),$
(54) $\displaystyle\frac{b_{\rm
R}(\theta)}{a+b}=\frac{1-r(\theta)}{2}\cos^{2}\left(\frac{\phi_{b}(\theta)}{2}\right).$
(55)

$\phi_{a}$ and $\phi_{b}$ are real and $r$ satisfies $-1\leq r(\theta)\leq 1$.
Since $a+b$ determines the overall time scale, we take it constant. Then, the
dynamical current is written as

$\displaystyle J_{\rm
d}=\frac{a+b}{8}\int_{0}^{2\pi}\frac{d\theta}{2\pi}\,(1-r^{2}(\theta))(\cos\phi_{a}(\theta)-\cos\phi_{b}(\theta)).$

It is nonzero when $\phi_{a}\neq\phi_{b}$ as we can understand from the
expression of the affinity at each time:

$\displaystyle{\cal A}(\theta)=\ln\frac{b_{\rm L}(\theta)a_{\rm
R}(\theta)}{a_{\rm L}(\theta)b_{\rm
R}(\theta)}=\ln\left(\frac{\tan\frac{\phi_{b}(\theta)}{2}}{\tan\frac{\phi_{a}(\theta)}{2}}\right)^{2}.$
(57)

On the other hand, the geometrical current is given by

$\displaystyle J_{\rm g}$ $\displaystyle=$
$\displaystyle\omega\int_{0}^{2\pi}\frac{d\theta}{2\pi}\,\frac{1}{8}\frac{dr(\theta)}{d\theta}$
$\displaystyle\times\left[(1+r(\theta))\cos\phi_{a}(\theta)+(1-r(\theta))\cos\phi_{b}(\theta)\right].$

To find a nonzero value of $J_{\rm g}$, we need that $r$ and, at least, one of
$\phi_{a}$ and $\phi_{b}$ has a nontrivial time dependence.

In a typical parametrization, the dynamical part is the dominant contribution
to the current. The dynamical current is calculated from the instantaneous
eigenvalue of the transition-rate matrix. The corresponding approximation in
the discrete-time case is given by

$\displaystyle
g(\chi)=\frac{1}{N}\ln(1+\Lambda_{0}^{\chi})\sim\frac{1}{N}\sum_{k=1}^{N}\ln(1+\Lambda_{0}^{\chi}(k)),$
(59)

where $\Lambda_{0}^{\chi}(k)$ denotes the largest eigenvalue of $W^{\chi}(k)$.
The geometrical current is typically a small quantity and becomes important
only when the dynamical current is absent.

Figure 4: Two protocols used in Sec. IV. $J_{\rm d}\neq 0$ and $J_{\rm g}=0$
for the protocol in the panel (a), and $J_{\rm d}=0$ and $J_{\rm g}\neq 0$ in
the panel (b).

The present parametrization allows us to study each part of the current
separately. In the first protocol, we use

$\displaystyle r(\theta)=\frac{1}{2},$ (60)
$\displaystyle\phi_{a}(\theta)=\frac{3\pi}{4}+\frac{\pi}{5}\cos\theta,$ (61)
$\displaystyle\phi_{b}(\theta)=\frac{\pi}{2}+\frac{2\pi}{5}\sin\theta.$ (62)

The corresponding behavior of parameters $a_{\rm L}$, $b_{\rm L}$, $a_{\rm
R}$, and $b_{\rm R}$ is shown in Fig. 4(a). In this case, $J_{\rm g}=0$ in the
continuous-time case and the dominant contribution of the current comes from
the dynamical part.

Figure 5: The density of zeros in the first protocol of Fig. 4(a) where
$J_{\rm d}\neq 0$ and $J_{\rm g}=0$ in the corresponding continuous-time
process. We set $N=2$ in the panel (a) and $N=4$ in the panel (b).  Figure 6:
$\varphi(J)$ in the first protocol. We set $N=2$ in the panel (a) and $N=4$ in
the panel (b). In each panel, the solid line denotes the result calculated
from Eq. (40). We use the density of zeros in Fig. 5. The dashed line denotes
the approximation in Eq. (59). The marked points denote simulation results. We
set $MN=5000$ and the number of samples to $10^{8}$ in the simulations.

In the discrete-time case, we consider the process by using the discretized
protocol. In Fig. 5, we plot the density of zeros at $N=2$ and $N=4$. We find
that each $z_{i}$ ($i=1,\dots,2N$) in these cases takes a negative value.
Then, the zeros continuously distribute on the negative real axis of $z$ with
$z_{i}\geq z\geq z_{i+1}$ ($i=1,2,\dots,2N-1$) where we set
$0>z_{1}>z_{2}>\cdots>z_{2N}$.

The corresponding rate function $\varphi(J)$, calculated from Eq. (40), is
shown in Fig. 6. We compare three results, $\varphi(J)$ exactly calculated
from $\rho(z)$, $\varphi(J)$ from the approximation in Eq. (59), and
$\varphi(J)$ from the simulations. As we see in Fig. 6, the result from zeros
is consistent with that from the simulations. We also observe that the average
current is well described by the approximation in Eq. (59).

We next consider the second protocol:

$\displaystyle r(\theta)=\frac{1}{2}+\frac{2}{5}\sin\theta,$ (63)
$\displaystyle\phi_{a}(\theta)=\phi_{b}(\theta)=\frac{\pi}{2}+\frac{\pi}{5}\cos\theta.$
(64)

The corresponding behavior of parameters $a_{\rm L}$, $b_{\rm L}$, $a_{\rm
R}$, and $b_{\rm R}$ is shown in Fig. 4(b). In this case, the current in the
continuous-time case is purely geometric. ${\cal A(\theta)}=0$ at each time,
which means that the instantaneous dynamical current is exactly equal to zero.
On the other hand, the flux in Eq. (LABEL:jg) is nonzero and we can observe a
nonzero geometrical current.

Figure 7:  The density of zeros in the second protocol of Fig. 4(b) where
$J_{\rm d}=0$ and $J_{\rm g}\neq 0$ in the corresponding continuous-time
process. We set $N=2$ in the panel (a) and $N=4$ in the panel (b).  Figure 8:
$\varphi(J)$ in the second protocol. We use the density of zeros in Fig. 7.
See the caption of Fig. 6 for other remarks.

We study how the geometric effect appears in the discrete-time case. At $N=2$
and $N=4$, each $z_{i}$ ($i=1,\dots,2N$) takes a negative value also in this
case. The results of $\rho(z)$ are shown in Fig. 7 and $\varphi(J)$ in Fig. 8.

In the case of $N=2$, the current distribution is almost symmetric and the
magnitude of the average current takes a small value. $\varphi(J)$ at small
$|J|$ is well described by the approximation in Eq. (59). At $N=2$, the
protocol trajectory in parameter space does not enclose any surface and no
geometric effect is observed. This result is significantly changed when we
consider the $N=4$ case. The approximation in Eq. (59) does not give any
reasonable result and we observe a nonzero average current. We find that
simulation results are well fitted by the result from zeros. Since the
condition of the adiabatic approximation is not obvious in the discrete-time
case, it is remarkable to find that the geometric effect can be seen at small
$N$.

##  V Floquet–Magnus expansion

We have described the systems with small frequency $\omega$ by using the
adiabatic approximation in the previous section. A similar systematic analysis
is possible when $\omega$ takes a large value. When the transition-rate matrix
drastically changes in short times, it is not appropriate to treat the
discrete-time system and we consider the continuous-time case in this section.

In contrast to the adiabatic case, we can easily find two zeros in two-state
systems and study the fluctuation effects on the current distributions. When
the period of the modulation is very small, we can expect that the system is
basically described by the averaged transition-rate matrix
$\bar{w}^{\chi}\sim\bar{w}_{0}^{\chi}=\int_{0}^{2\pi}\frac{d\theta}{2\pi}w^{\chi}(\theta)$.
The correction can be expressed by using a series expansion of the inverse
frequency, which is known as the Floquet–Magnus expansion [23, 24].

We define the effective transition-rate matrix $\bar{w}^{\chi}$ from the
relation

$\displaystyle\exp\left(\frac{2\pi}{\omega}\bar{w}^{\chi}\right)={\rm
T}\exp\left(\int_{0}^{2\pi/\omega}dt\,w^{\chi}(\omega t)\right),$ (65)

where T denotes that the time-evolution operator is represented by the time-
ordered product. In the Floquet–Magnus expansion, this effective matrix is
represented as $\bar{w}^{\chi}=\sum_{k=0}^{\infty}\bar{w}_{k}^{\chi}$.
$\bar{w}_{k}^{\chi}$ is proportional to $1/\omega^{k}$ and the explicit forms
at $k=1$ and 2 are given as follows:

$\displaystyle\bar{w}_{1}^{\chi}=\frac{1}{2\omega}\int_{0}^{2\pi}\frac{d\theta_{1}}{2\pi}\int_{0}^{\theta_{1}}d\theta_{2}[w^{\chi}(\theta_{1}),w^{\chi}(\theta_{2})],$
(66) $\displaystyle\bar{w}_{2}^{\chi}$ $\displaystyle=$
$\displaystyle\frac{1}{6\omega^{2}}\int_{0}^{2\pi}\frac{d\theta_{1}}{2\pi}\int_{0}^{\theta_{1}}d\theta_{2}\int_{0}^{\theta_{2}}d\theta_{3}$
(67)
$\displaystyle([w^{\chi}(\theta_{1}),[w^{\chi}(\theta_{2}),w^{\chi}(\theta_{3})]$
$\displaystyle+[w^{\chi}(\theta_{3}),[w^{\chi}(\theta_{2}),w^{\chi}(\theta_{1})]).$

To give a concrete discussion, we treat the two-state case. Since the
expansion keeps the trace of the transition-rate matrix invariant, we can
generally write

$\displaystyle\bar{w}^{\chi}=-\frac{\bar{a}+\bar{b}}{2}+\left(\begin{array}[]{cc}\alpha^{\chi}&\beta^{\chi}\\\
\gamma^{\chi}&-\alpha^{\chi}\end{array}\right),$ (70)

where

$\displaystyle\bar{a}=\int_{0}^{2\pi}\frac{d\theta}{2\pi}a(\theta),$ (71)

and $\bar{b}$ is defined in a similar way. The cumulant generating function is
given by the largest eigenvalue of $\bar{w}^{\chi}$ as

$\displaystyle
g(\chi)=-\frac{\bar{a}+\bar{b}}{2}+\sqrt{(\alpha^{\chi})^{2}+\beta^{\chi}\gamma^{\chi}}.$
(72)

Up to the first order of the expansion, each element is given respectively by

$\displaystyle\alpha^{\chi}$ $\displaystyle\sim$
$\displaystyle\frac{\bar{a}-\bar{b}}{2}+\frac{1}{2\omega}\int_{0}^{2\pi}\frac{d\theta_{1}}{2\pi}\int_{0}^{\theta_{1}}d\theta_{2}\left(a^{\chi}(\theta_{1})b^{\chi}(\theta_{2})\right.$
(73) $\displaystyle\left.-b^{\chi}(\theta_{1})a^{\chi}(\theta_{2})\right),$
$\displaystyle\beta^{\chi}$ $\displaystyle\sim$
$\displaystyle\bar{a}^{\chi}+\frac{1}{2\omega}\int_{0}^{2\pi}\frac{d\theta_{1}}{2\pi}\int_{0}^{\theta_{1}}d\theta_{2}[(a(\theta_{1})-b(\theta_{1}))a^{\chi}(\theta_{2})$
(74) $\displaystyle-a^{\chi}(\theta_{1})(a(\theta_{2})-b(\theta_{2}))],$
$\displaystyle\gamma^{\chi}$ $\displaystyle\sim$
$\displaystyle\bar{b}^{\chi}+\frac{1}{2\omega}\int_{0}^{2\pi}\frac{d\theta_{1}}{2\pi}\int_{0}^{\theta_{1}}d\theta_{2}[-(a(\theta_{1})-b(\theta_{1}))b^{\chi}(\theta_{2})$
(75) $\displaystyle+b^{\chi}(\theta_{1})(a(\theta_{2})-b(\theta_{2}))].$

From these relations, we can find that the cumulant generating function takes
a form

$\displaystyle
g(\chi)\sim-\frac{\bar{a}+\bar{b}}{2}+\frac{\bar{a}+\bar{b}}{2}\sqrt{\frac{(e^{\chi}-z_{1})(e^{\chi}-z_{2})}{e^{\chi}(1-z_{1})(1-z_{2})}}.$
(76)

As in the zeroth order, the function is characterized by two zeros.

We can also examine the higher-order corrections. Although it is a difficult
task to find the explicit forms, the expansion is systematically represented
by multiple commutators. We find that up to the $2n$th order of the expansion

$\displaystyle
g(\chi)\sim-\frac{\bar{a}+\bar{b}}{2}+\frac{\bar{a}+\bar{b}}{2}\sqrt{\prod_{k=0}^{n-1}\frac{(e^{\chi}-z_{1}(k))(e^{\chi}-z_{2}(k))}{e^{\chi}(1-z_{1}(k))(1-z_{2}(k))}}.$

Then, the average current is given by

$\displaystyle
J\sim\frac{\bar{a}+\bar{b}}{2}\sum_{k=0}^{n-1}\frac{1}{2}\left(\frac{1}{1-z_{1}(k)}+\frac{1}{1-z_{2}(k)}-1\right).$
(78)

The current is generally characterized by many zeros and the number of the
zeros increases as we increase the order of the expansion.

Figure 9: Protocol used in Sec. V. (a) $a_{\rm L}$, $a_{\rm R}$, and $a=a_{\rm
L}+a_{\rm R}$. (b) $b_{\rm L}$, $b_{\rm R}$, and $b=b_{\rm L}+b_{\rm R}$.
Figure 10: Zeros obtained by the Floquet–Magnus expansion. At the first-order
expansion, we have two zeros at each frequency $\omega$, which are denoted by
points $\circ$. At the second order, we have four zeros denoted by points
$\times$ and one of them is shown in the inset.  Figure 11: Average current as
a function of the inverse frequency. The bold solid line denotes the exact
result, dotted line denotes the first-order accuracy of the Floquet–Magnus
expansion, dashed line denotes the second order of the expansion.

We show the result of the expansion up to the second order. We set the
protocol in Eqs. (52)-(55) with

$\displaystyle r=\frac{\pi}{4}\sin\theta,$ (79)
$\displaystyle\phi_{a}=\frac{\pi}{2}+\frac{\pi}{4}\cos\theta,$ (80)
$\displaystyle\phi_{b}=\frac{\pi}{2}+\frac{\pi}{4}\sin\theta,$ (81)

as shown in Fig. 9. The zeros are plotted in Fig. 10 and the corresponding
current is in Fig. 11. We see that the current is properly described by using
the zeros. At the zeroth order of the expansion, two zeros are located on the
negative real axis. They change as functions of $1/\omega$ by taking into
account the higher-order corrections. At the same time, when we increase the
order of the expansion, different zeros appear around the origin and on the
real axis far from the origin. These zeros are treated in pairs and give small
corrections to the current as

$\displaystyle\frac{1}{1-z_{1}(k)}+\frac{1}{1-z_{2}(k)}-1\to 0,$ (82)

at $|z_{1}(k)|\to\infty$ and $z_{2}(k)\to 0$.

To study the meaning of zeros, we examine how the fluctuation theorem is
affected by the expansion. The rate is calculated from the Legendre
transformation $\varphi(J)=g(\chi)-J\chi$. In the expansion up to $2n$th
order, the relation between $\chi$ and $J$ is obtained from

$\displaystyle J(\chi)$ $\displaystyle\sim$
$\displaystyle\frac{\bar{a}+\bar{b}}{2}\sum_{k=0}^{n-1}\frac{1}{2}\left(\frac{1}{e^{\chi}-z_{1}(k)}+\frac{1}{e^{\chi}-z_{2}(k)}-\frac{1}{e^{\chi}}\right)$
(83)
$\displaystyle\times\sqrt{\prod_{k=0}^{n-1}\frac{(e^{\chi}-z_{1}(k))(e^{\chi}-z_{2}(k))}{e^{\chi}(1-z_{1}(k))(1-z_{2}(k))}}.$

$\chi$ with $e^{\chi}>0$ is uniquely determined for a given $J$. Since it is
difficult to obtain a compact form of the general solution, we discuss the
cases where $J$ takes extremal values.

When the absolute value of $J$ is large, $|J|\gg\bar{a}+\bar{b}$, the
corresponding counting field is also large $|\chi|\gg 1$ and we can evaluate

$\displaystyle
J(\chi)\sim\frac{\bar{a}+\bar{b}}{2}n(e^{\chi})^{n/2}\sqrt{\prod_{k=0}^{n-1}\frac{1}{(1-z_{1}(k))(1-z_{2}(k))}}$

for $e^{\chi}\gg 1$ and

$\displaystyle
J(\chi)\sim-\frac{\bar{a}+\bar{b}}{2}n(e^{\chi})^{-n/2}\sqrt{\prod_{k=0}^{n-1}\frac{z_{1}(k)z_{2}(k)}{(1-z_{1}(k))(1-z_{2}(k))}}$

for $e^{\chi}\ll 1$. Then, after some calculations we obtain

$\displaystyle\varphi(J)-\varphi(-J)\sim
J\ln\left(\prod_{k=0}^{n-1}z_{1}(k)z_{2}(k)\right)^{1/n}$ (86)

for $|J|\gg\bar{a}+\bar{b}$. This relation implies that the affinity is given
by the geometric mean of zeros. As we find in the above example, the product
$z_{1}(k)z_{2}(k)$ takes a finite value.

In the opposite limit at $J=0$, the corresponding $e^{\chi}$ is obtained by
solving

$\displaystyle\sum_{k=0}^{n-1}\frac{(e^{\chi^{*}})^{2}-z_{1}(k)z_{2}(k)}{(e^{\chi^{*}}-z_{1}(k))(e^{\chi^{*}}-z_{2}(k))}=0.$
(87)

Then, $\varphi(J)$ is obtained as $\varphi(J)\sim J\ln e^{\chi^{*}}$ and

$\displaystyle\varphi(J)-\varphi(-J)\sim J\ln(e^{\chi^{*}})^{2},$ (88)

for $|J|\ll\bar{a}+\bar{b}$.

We numerically find that the product $z_{1}(k)z_{2}(k)$ does not satisfy Eq.
(35) and is dependent on the index $k$, which means that Eq. (86) does not
hold in the time-dependent system. The fluctuation theorem still holds in that
case, but its consequence is not represented in a simple form as in Eq. (86)
[26].

##  VI Summary

We have studied dynamical Lee–Yang zeros in stochastic processes. As we
mentioned in the Introduction, the distribution of zeros has complete
information about the statistical properties of the processes. By discretizing
the time variable, we can explicitly confirm that the number of zeros is
related not to the number of states but to the number of events. We find
nontrivial distributions of zeros even in simple two-state processes. The
distribution of zeros is related to the current distributions and we expect
that our discrete-time result can be confirmed experimentally by measuring
processes with rare events.

Since our analysis is restricted to simple two-state systems, the zeros
distribute only on the negative real axis of $z=e^{\chi}$. In statistical
mechanics, thermal phase transitions can be found by studying the distribution
of zeros. The locations of zeros represent the boundary of different phases in
complex plane. We have not found closed boundaries in the present analysis. It
may be an interesting problem to study more complex systems such as quantum
dot systems coupled with external reservoirs treated in previous studies [14,
18].

One of the important findings in the present study is that the method can be
applied in the case of the periodically oscillating transition rate. In the
adiabatic case with slowly varying transition-rate matrix, the property of the
density of zeros is dependent on the geometric properties of the protocols. We
also find in the fast-driving regime of the continuous-time case that the
number of zeros depends on the order of the expansion. These results show that
it is important to study global distributions of zeros to characterize the
dynamically fluctuating systems instead of using a few number of zeros. The
distribution also reflects a fundamental symmetry of the system. We expect
that the principle of nonequilibrium thermodynamics can be studied along with
the method developed in the present study.

## Acknowledgments

We thank Y. Utsumi for valuable comments. K.T. was supported by JSPS KAKENHI
Grants No. JP20K03781 and No. JP20H01827.

## References

## References

  * [1] D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Probability of Second Law Violations in Shearing Steady States,  Phys. Rev. Lett. 71, 2401 (1993).
  * [2] G. Gallavotti and E. G. D. Cohen, Dynamical Ensembles in Nonequilibrium Statistical Mechanics,  Phys. Rev. Lett. 74, 2694 (1995).
  * [3] G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences,  Phys. Rev. E 60, 2721 (1999).
  * [4] L. S. Levitov and G. B. Lesovik, Charge distribution in quantum shot noise, JETP Lett. 58, 230 (1993).
  * [5] L. S. Levitov, H.W. Lee, and G. B. Lesovik, Electron counting statistics and coherent states of electric current,  J. Math. Phys. 37, 4845 (1996).
  * [6] D. A. Bagrets and Y. V. Nazarov, Full counting statistics of charge transfer in Coulomb blockade systems,  Phys. Rev. B 67, 085316 (2003).
  * [7] K. Saito and Y. Utsumi, Symmetry in full counting statistics, fluctuation theorem, and relations among nonlinear transport coefficients in the presence of a magnetic field,  Phys. Rev. B 78, 115429 (2008).
  * [8] C. N. Yang and T. D. Lee, Statistical theory of equations of state and phase transitions. I. Theory of condensation,  Phys. Rev. 87, 404 (1952).
  * [9] T. D. Lee and C. N. Yang, Statistical theory of equations of state and phase transitions. II. Lattice gas and Ising model,  Phys. Rev. 87, 410 (1952).
  * [10] M. E. Fisher, The Nature of Critical Points, Lectures in Theoretical Physics Vol. VII C (University of Colorado Press, Boulder, 1965).
  * [11] A. G. Abanov and D. A. Ivanov, Allowed Charge Transfers between Coherent Conductors Driven by a Time-Dependent Scatterer,  Phys. Rev. Lett. 100, 086602 (2008).
  * [12] D. A. Ivanov and A. G. Abanov, Phase transitions in full counting statistics for periodic pumping,  Europhys. Lett. 92, 37008 (2010).
  * [13] C. Flindt and J. P. Garrahan, Trajectory Phase Transitions, Lee–Yang Zeros, and High-Order Cumulants in Full Counting Statistics,  Phys. Rev. Lett. 110, 050601 (2013).
  * [14] Y. Utsumi, O. Entin-Wohlman, A. Ueda, and A. Aharony, Full-counting statistics for molecular junctions: Fluctuation theorem and singularities,  Phys. Rev. B 87, 115407 (2013).
  * [15] J. M. Hickey, C. Flindt, and J. P. Garrahan, Intermittency and dynamical Lee–Yang zeros of open quantum systems,  Phys. Rev. E 90, 062128 (2014).
  * [16] R. S. Souto, A. Martín-Rodero, and A. L. Yeyati, Analysis of universality in transient dynamics of coherent electronic transport,  Fortschr. Phys. 65, 1600062 (2016).
  * [17] K. Brandner, V. F. Maisi, J. P. Pekola, J. P. Garrahan, and C. Flindt, Experimental Determination of Dynamical Lee–Yang Zeros,  Phys. Rev. Lett. 118, 180601 (2017).
  * [18] R. S. Souto, A. Martín-Rodero, and A. L. Yeyati, Quench dynamics in superconducting nanojunctions: Metastability and dynamical Yang–Lee zeros,  Phys. Rev. B 96, 165444 (2017).
  * [19] B.-B. Wei and R.-B. Liu, Lee–Yang Zeros and Critical Times in Decoherence of a Probe Spin Coupled to a Bath,  Phys. Rev. Lett. 109, 185701 (2012).
  * [20] X. Peng, H. Zhou, B.-B. Wei, J. Cui, J. Du, and R.-B. Liu, Experimental Observation of Lee–Yang Zeros,  Phys. Rev. Lett. 114, 010601 (2015).
  * [21] D. J. Thouless, Quantization of particle transport,  Phys. Rev. B 27, 6083 (1983).
  * [22] N. A. Sinitsyn and I. Nemenman, The Berry phase and the pump flux in stochastic chemical kinetics,  Europhys. Lett. 77, 58001 (2007).
  * [23] W. Magnus, On the exponential solution of differential equations for a linear operator,  Commun. Pure Appl. Math. 7, 649 (1954).
  * [24] S. Blanes, F. Casas, J. Oteo, and J. Ros, The Magnus expansion and some of its applications,  Phys. Rep. 470, 151 (2009).
  * [25] K. Takahashi, K. Fujii, Y. Hino, and H. Hayakawa, Nonadiabatic Control of Geometric Pumping,  Phys. Rev. Lett. 124, 150602 (2020).
  * [26] K. Takahashi, Y. Hino, K. Fujii, and H. Hayakawa, Full counting statistics and fluctuation-dissipation relation for periodically driven two-state systems,  J. Stat. Phys. 181, 2206 (2020).
  * [27] M. V. Berry, Quantal phase factors accompanying adiabatic changes,  Proc. R. Soc. London, Ser. A 392, 45 (1984).


148 2021 submission date publication date NeurIPS 2020 Preregistration
Workshop

# Testing the Genomic Bottleneck Hypothesis in Hebbian Meta-Learning

Rasmus Berg Palm<EMAIL_ADDRESS>
Elias Najarro<EMAIL_ADDRESS>
Sebastian Risi<EMAIL_ADDRESS>
IT University of Copenhagen

###### Abstract

Hebbian meta-learning has recently shown promise to solve hard reinforcement
learning problems, allowing agents to adapt to some degree to changes in the
environment. However, because each synapse in these approaches can learn a
very specific learning rule, the ability to generalize to very different
situations is likely reduced. We hypothesize that limiting the number of
Hebbian learning rules through a “genomic bottleneck” can act as a regularizer
leading to better generalization across changes to the environment. We test
this hypothesis by decoupling the number of Hebbian learning rules from the
number of synapses and systematically varying the number of Hebbian learning
rules. The results in this paper suggest that simultaneously learning the
Hebbian learning rules and their assignment to synapses is a difficult
optimization problem, leading to poor performance in the environments tested.
However, parallel research to ours finds that it is indeed possible to reduce
the number of learning rules by clustering similar rules together. How to best
implement a “genomic bottleneck” algorithm is thus an important research
direction that warrants further investigation.

###### keywords:

Hebbian learning, Meta learning, Genomic Bottleneck

##  1 Introduction

Deep reinforcement learning has made great progress recently on very hard
problems like Go, Starcraft and Dota using deep neural networks (Silver et
al., 2016; Vinyals et al., 2019; Berner et al., 2019). However, once learned,
the networks are static and highly specific. As such, there is very little
capacity to adapt to changes in the environment or to generalize across
environments (Justesen et al., 2018). For instance, the state-of-the-art
AlphaStar agent trained to play one race in Starcraft cannot play another,
even though it is a very similar task, and much less play Go or load a
dishwasher. In contrast animals and humans show remarkable flexibility in
their ability to generalize across tasks and adapt to changes.

Meta-learning proposes to overcome these limitations by learning-to-learn.
That is to learn general learning rules that are broadly applicable and enable
an agent to quickly adapt to changes in the environment or new tasks (Finn et
al., 2017; Wang et al., 2016; Weng, 2018).

One particularly interesting approach to meta-learning is Hebbian meta-
learning. The goal in Hebbian meta-learning is to learn Hebbian learning rules
that enable an agent to quickly learn to perform well in its environment and
adapt to changes. Hebbian learning is a learning paradigm in which the
synaptic strength between neurons is determined by the correlation of their
activity (Hebb, 2005); informally, “neurons that fire together, wire
together”. The Hebbian learning paradigm is promising since it proposes a
simple and _general_ learning paradigm supported by extensive empirical
evidence in biological brains (Prezioso et al., 2018; Caporale and Dan, 2008).

Najarro and Risi (2020) showed promising results using Hebbian meta-learning
to solve a difficult car racing reinforcement problem and quadruped
locomotion. Impressively the agent’s policy network was initialized randomly
each episode, so it had to learn the task during its lifetime using only the
meta-learned Hebbian learning rules. Further, the learned Hebbian learning
rules generalized to an unseen quadruped leg damage scenario, which baseline
non-plastic feedforward neural networks could not. However, each synapse had
its own learning rule, which allowed the network to learn very specific
learning rules for each neuron. This raises the question to which degree the
network learned to learn, or whether each learning rule rather encoded very
specific dynamics for each synapse.

We hypothesize that the architecture in Najarro and Risi (2020) is limited in
its ability to discover general learning rules, but rather encodes very
specific dynamics for each synapse. We further hypothesize that in order to
learn generally useful learning rules the network must be limited to how many
learning rules it can learn. This is inspired by the genomic bottleneck
observed in humans and complex animals, where the information that can be
stored in the genome is several orders of magnitude smaller than what is
needed to determine the final wiring of the brain (Zador, 2019). We
hypothesize that limiting the number of learning rules acts as a regularizer,
which improves generalization and adaptability.

##  2 Related Work

Meta-Learning. In meta-learning, the goal is to learn to quickly adapt to a
target task given a set of training tasks (Weng, 2018). Common approaches can
loosely be categorized into black-box, optimization, and metric-based. Black-
box methods learn a function that, conditioned on samples from a new task,
outputs a function to solve the new task. For instance by jointly learning a
network that can produce a latent summary of a new task and a network that can
solve the task given the latent summary (Santoro et al., 2016; Mishra et al.,
2017). Optimization-based attempts to learn to quickly optimize on a new task,
e.g. by finding an initialization from which optimization on new tasks is fast
(Finn et al., 2017; Nichol et al., 2018) or by learning the optimization
algorithm (Andrychowicz et al., 2016; Li et al., 2017). Metric based
approaches learns an embedding space that facilitates effective distance-based
classification, e.g. Siamese networks (Koch et al., 2015) or prototypical
networks (Snell et al., 2017). Meta Reinforcement Learning (meta-RL) extends
the meta-learning idea to the reinforcement learning setting. Formally the
goal is to quickly learn to perform well in a new Markov Decision Process
(MDP) given a set of training MDPs. A common approach is to learn a recurrent
neural network policy where the hidden activations are not reset between
episodes (Wang et al., 2016; Duan et al., 2016). This allows the policy
network to discover how the environment behaves across episodes and adapt its
policy. Another common approach is meta-learning an initialization of a policy
network from which policy gradient descent can quickly adapt to the new MDP
(Finn et al., 2017; Song et al., 2019).

Plastic Artificial Neural Networks A less explored meta-learning approach is
based on plastic neural networks that are optimized to have both innate
properties and the ability to learn during their lifetime. For example, such
networks can learn by changing the connectivity among neurons through local
learning rules like Hebbian plasticity (Soltoggio et al., 2018, 2007). Often
these networks are optimized through evolutionary algorithms (Soltoggio et
al., 2018; Najarro and Risi, 2020) but more recently optimizing the plasticity
of connections in a network through gradient descent has also been shown
possible (Miconi et al., 2018). However, in contrast to the evolving Hebbian
learning rules approach in Najarro and Risi (2020), the gradient descent
approach (Miconi et al., 2018) was so far restricted to only evolving a single
plasticity parameter for each connection instead of a different Hebbian rule.

Fast Weights Artificial neural networks (ANN) have either a slow form of
storing information—through updating the weights—or, if they are recurrent, a
very fast form of information storage in the form of internal activations.
Fast weights seek to introduce an intermediate time-scale to the information
storage in ANNs (Hinton and Plaut, 1987; Schmidhuber, 1992) and is motivated
by the observation that biological neural networks have learning processes
occurring concurrently that span across very different time scales. Recently,
this approach has been successfully applied to image recognition tasks as well
as a model for content-addressable memory (Ba et al., 2016). In the context of
meta-learning, fast weights has been shown to perform meta-RL by having an ANN
update the weights of a policy network on a per-task basis (Munkhdalai and Yu,
2017). Additionally, Munkhdalai and Trischler (2018) showed that a fast-
weights Hebbian mechanism is capable of performing one-shot supervised
learning tasks.

##  3 Hebbian meta-learning

In Hebbian meta-learning, the goal is to meta-learn Hebbian learning rules
that enable an agent to perform well, and adapt to changes, in its
environment.

###  3.1 Hebbian Learning

The agent acts in an episodic reinforcement learning environment and is
controlled by a neural policy network. At the start of each episode, the
policy network is initialized with random weights, which then undergoes
changes according to Hebbian learning rules during the episode. Specifically,
we use the ABCD Hebbian weight plasticity formalization (Soltoggio et al.,
2007) and the weights are updated at each step of the episode. The change to
the weight connecting neuron $i$ and $j$ is

$\displaystyle\Delta
w_{ij}=\eta_{ij}\left(A_{ij}o_{i}o_{j}+B_{ij}o_{i}+C_{ij}o_{j}+D_{ij}\right)\,,$
(1)

where $o_{i}$ and $o_{j}$ are the pre- and post-synaptic activation of the
neurons and $h=\\{\eta,A,B,C,D\\}$ are the Hebbian parameters that are fixed
during the episode.

###  3.2 Evolutionary Strategies

The Hebbian parameters are meta-learned on an evolutionary time scale $t$,
using evolutionary strategies (ES) (Salimans et al., 2017). In evolutionary
strategies the goal is to maximize the expected fitness of a distribution of
individuals, $\max_{\theta}\mathbb{E}_{z\sim p(z|\theta)}F(z)$, where $F(z)$
is the fitness of an individual $z$ and $\theta$ parameterize this
distribution. We compute the gradient of this objective using the score
function estimator and use gradient ascent to maximize it,

$\displaystyle\nabla_{\theta}\mathbb{E}_{z\sim p(z|\theta)}F(z)$
$\displaystyle=\mathbb{E}_{z\sim
p(z|\theta)}F(z)\nabla_{\theta}\log(p(z|\theta))$ $\displaystyle\theta^{t+1}$
$\displaystyle=\theta^{t}+\alpha\left[\mathbb{E}_{z\sim
p(z|\theta^{t})}F(z)\nabla_{\theta^{t}}\log(p(z|\theta^{t}))\right]\,,$ (2)

where $\alpha$ is the learning rate. The expectation is evaluated using $n$
samples, the population size. In order to use ES then, we must define $F(z)$
and $p(z|\theta)$. In this paper $F(z)$ is always the accumulated reward of an
episode.

###  3.3 Individual Learning Rules

For the case of individual learning rules, $N$ synapses each have their own
learning rules $h_{i}=[\eta_{i},A_{i},B_{i},C_{i},D_{i}]$, which are drawn
from independent normal distributions with fixed variance $\sigma^{2}$ and
meta-learned means $\mu\in\mathbb{R}^{N\times 5}$. In this case
$p(z|\theta)=p(h|\mu)=\prod_{i=1}^{N}\mathcal{N}(h_{i}|\mu_{i},\sigma)$, where
$\mu_{i}$ denotes the $i$’th row of the $\mu$ parameters. Inserting into eq.
(2) and deriving the gradient this reduces to the expression in (Najarro and
Risi, 2020; Salimans et al., 2017),

$\displaystyle\mu^{t+1}=\mu^{t}+\alpha\left[\frac{1}{\sigma}\mathbb{E}_{\epsilon\sim\mathcal{N}(0,I)}F(\mu^{t}+\sigma\epsilon)\epsilon\right]\,,$

where $\epsilon$ is a sample from a standard normal $\mathcal{N}(0,I)$.

###  3.4 Shared Learning Rules

To explore the effect of sharing a limited amount of Hebbian learning rules we
sample each synapses’ learning rule from a Gaussian Mixture Model (GMM) with
$M$ clusters, such that

$\displaystyle k$ $\displaystyle\sim
p(k|\lambda)\in[1,...,M]^{N}\,,\lambda\in\mathbb{R}^{N\times M}\,,$
$\displaystyle h$
$\displaystyle\sim\mathcal{N}(h|\mu_{k},\sigma)\in\mathbb{R}^{N\times
5}\,,\mu\in\mathbb{R}^{M\times 5}\,,$ $\displaystyle
p(z|\theta)=p(h|\mu,\lambda)$
$\displaystyle=\prod_{i=1}^{N}\sum_{k=1}^{M}\mathcal{N}(h_{i}|\mu_{k},\sigma)p(k|\lambda_{i})\,.$
(3)

Here $p(k|\lambda_{i})=e^{\lambda_{ik}}/\sum_{j=1}^{M}e^{\lambda_{ij}}$, i.e.
the softmax categorical distribution parameterized by $\lambda_{i}$ logits for
synapse $i$ and the meta learned parameters are $\mu$ and $\lambda$. We use
automatic differentiation to compute the gradient of the log likelihood of
this distribution with respect to $\mu$ and $\lambda$, and then use eq. (2) to
update $\mu$ and $\lambda$.

This approach to assigning learning rules to synapses is the most flexible and
direct but requires $N\times M$ parameters. However, we are only interested in
testing the effect of limiting the number of learning rules, not the number of
parameters or bits needed to encode an individual, although those are
interesting directions for future work.

##  4 Experiments

The first experiment is to replicate the original results of Najarro and Risi
(2020) to ensure a fair comparison. In all experiments, we use the same
experimental setup, architectures, and hyper-parameters as in Najarro and Risi
(2020) except where noted.

Similar to Najarro and Risi (2020) we experiment on the car racing and
quadruped locomotion tasks. We perform leave-one-out cross-validation, with
five variations for each task. The five variations for the car racing tasks
are (1) default settings, (2,3) twice and half the road friction coefficient,
and (4,5) constant force pushing the car west and east. For the quadruped
locomotion task, the variations are (1) default settings, (2,3) left and right
front leg damage as in Najarro and Risi (2020) and (4,5) 50% longer rear and
front legs. We perform leave-one-out cross-validation by leaving one variation
out for testing and training on the remaining variations. Specifically, we use
eq. (2) to maximize $\mathbb{E}_{z\sim p(z|\theta)}F^{\prime}(z)$ where
$F^{\prime}$ is the fitness function of a meta-task formed by taking the
average fitness across the four training variations of the task.

To test our core hypothesis we evaluate for a varying number of learning rules
expressed as a fraction of the number of synapses such that
$M=\frac{1}{\rho}N$. We vary $\rho=[1,16,32,64,128,256,N]$, where $\rho=1$
corresponds to a learning rule per synapse and $\rho=N$ corresponds to a
single learning rule.

We compare to three baselines in all experiments: (1) A Hebbian meta-learning
network with a learning rule per synapse as in Najarro and Risi (2020), (2) an
identical static network with learned weights, and (3) a LSTM baseline with
hidden states initialized to zero at the start of each episode (Hochreiter and
Schmidhuber, 1997). We construct the LSTM baseline by replacing the first
densely connected layer in the static baseline with a LSTM layer and keeping
everything else the same. All architectures are optimized using ES.

##  5 Results

We use the EvoStrat computational library (Palm, 2020) for all experiments.
Code to reproduce all results are available at
github.com/rasmusbergpalm/hebbian-evolution.

###  5.1 Replication

We successfully replicated the results of Najarro and Risi (2020) for the car
racing and quadruped walker environments consistently achieving high rewards
in 300 and 500 generations respectively with a Hebbian network with a learning
rule per synapse.

###  5.2 Computational issues

Evaluating eq. 2 for the shared Hebbian learning rules for small values of
$\rho$, becomes very expensive in terms of memory due to the number of
parameters of the model. Unfortunately, that means we’re not able to perform
the experiments for values of $\rho<32$.

###  5.3 Ant Environment

The results on the quadruped walker environment "AntBulletEnv-v0" are
summarized in table 1. The networks are trained on the 4 left-most
morphologies and tested on the "Long front" morphology corresponding to 50%
longer front legs.
Baselines The baselines all do reasonably well. The static network is best on
average on the training tasks, whereas the LSTM and Hebbian learning rules are
better on the unseen test task, although the standard deviations are large.
This is in line with the original results from Najarro and Risi (2020).
Shared learning rules We first test with $\rho=128$, which fails on both the
training and testing morphology. Given that it fails on the training tasks we
only evaluate subsequent values of $\rho$ on the default training morphology.
We assume that if it fails here it will also fail under the more difficult
leave-one-out meta-learning testing scheme.

Table 1: Ant environment results. Average distance traveled. Bold indicate unseen testing morphology. Mean and standard deviation from 100 samples. Rounded to nearest integer. Model | Default | Damage left | Damage right | Long back | Long front
---|---|---|---|---|---
Static | $1545\pm 30$ | $839\pm 60$ | $1891\pm 105$ | $1197\pm 31$ | $202\pm 174$
LSTM | $1097\pm 57$ | $1374\pm 74$ | $1240\pm 140$ | $806\pm 113$ | $314\pm 115$
Hebbian | $928\pm 88$ | $1102\pm 164$ | $981\pm 140$ | $741\pm 169$ | $364\pm 131$
Shared($\rho=32$) | $2\pm 12$ |  |  |  |
Shared($\rho=64$) | $3\pm 13$ |  |  |  |
Shared($\rho=128$) | $0\pm 10$ | $-9\pm 18$ | $15\pm 31$ | $1\pm 6$ | $0\pm 4$
Shared($\rho=256$) | $2\pm 12$ |  |  |  |
Shared($\rho=N$) | $-1\pm 0$ |  |  |  |

###  5.4 Car Racing Environment

Baselines All the baselines learn to perform well, surprisingly also on the
unseen test environment variation.
Shared learning rules Similarly to the ant results, we test with $\rho=128$
first, and it fails to learn anything on both the training and testing tasks.
Given similar negative results from the ant environment, we conclude that
there is something fundamentally problematic about the shared learning rule
experiments, and focus on diagnosing that. As such we do not test for other
values of $\rho$ on the car environment.

Table 2: Car racing environment results. Average distance traveled. Model | Default | Half friction | Double friction | Pushed left | Pushed right
---|---|---|---|---|---
Static | $784\pm 119$ | $491\pm 202$ | $789\pm 202$ | $738\pm 219$ | $628\pm 228$
LSTM | $680\pm 247$ | $375\pm 209$ | $711\pm 243$ | $653\pm 273$ | $610\pm 316$
Hebbian | $674\pm 135$ | $552\pm 190$ | $692\pm 120$ | $725\pm 139$ | $606\pm 186$
Shared($\rho=128$) | $6\pm 3$ | $5\pm 1$ | $8\pm 7$ | $6\pm 1$ | $6\pm 2$

###  5.5 Post-hoc analysis

We deviate from our pre-planned experimental setup to see if we can get the
shared Hebbian learning rules to work. We experiment on the default Ant
morphology using $\rho=128$. We independently try lowering the learning rate
to $0.01$ (from $0.2$), reducing the number of shared learning rules to $16$
and initializing the shared learning rule means from $\mathcal{N}(0,10)$ (from
$\mathcal{N}(0,1)$). However, none of these modifications helped. We also try
randomly assigning shared learning rules instead of learning the assignments.
This modification improves performance, however, the networks are still much
worse than the baselines (not shown). We don’t pursue the random assignment
any further since we are interested in learning the assignments.

We note that the Hebbian networks with shared learning rules fail to learn the
mixing weights $\lambda$, which stay close to a uniform distribution, and that
the gradient estimates on the mixing weights are two orders of magnitude
smaller than on the learning rule means $\mu$. We manually verify that the
$\log p(h|\mu,\lambda,\sigma)$ computations and the gradients
$\nabla_{\mu}\log p(h|\mu,\lambda,\sigma)$ and $\nabla_{\lambda}\log
p(h|\mu,\lambda,\sigma)$ are correct for a simple test case.

Suspecting a difficult optimization landscape we use the Adam optimizer which
scales the learning rate individually for each parameter (Kingma and Ba,
2014). With the Adam optimizer, the networks attain poor but non-zero
performance. This rules out the most severe bugs in the implementation and
experimental setup and lends support to our hypothesis that it is “just” a
difficult optimization problem. We experiment with several learning rates and
learning rate schedules, achieving the best results using a learning rate of
$1.0$ and an exponential learning rate decay of $0.9931$ which corresponds to
halving the learning rate every 100 updates. However, these results are still
much worse than the baselines, achieving $555\pm 89$ reward on the default Ant
morphology.

##  6 Discussion

The experimental results in this paper do not support our initial hypothesis.
However, parallel research to ours found that it is indeed possible to reduce
the number of Hebbian learning rules significantly and improve generalization
on unseen environments by iteratively clustering similar learning rules and
re-training (Pedersen and Risi, 2021). As such, the results suggest that our
approach (i.e. modeling the learning rules as drawn from a Gaussian Mixture
Model and jointly optimizing the assignment and the learning rules) is a hard
optimization problem.

There is a bit of a chicken and egg problem in jointly optimizing the
assignment and the learning rules: it is challenging to optimize the learning
rules before they are assigned to synapses, and it is challenging to optimize
the assignment to synapses before the learning rules are learned. The approach
introduced by Pedersen and Risi (2021) avoids this issue by learning a Hebbian
learning rule per synapse, while periodically clustering similar learning
rules and re-training with the reduced set of learning rules. It could also be
possible to avoid this chicken and egg problem by alternately learning the
assignment and the learning rules in an Expectation-Maximization (EM) inspired
approach (Dempster et al., 1977).

Another possible explanation for the difficulty in directly learning a small
number of shared Hebbian learning rules is the lottery ticket hypothesis
(Frankle and Carbin, 2018). The lottery ticket hypothesis is based on the
observation that neural networks can often be pruned aggressively without
losing performance, and that it is in fact a small sub-network of the bigger
network which solves the task; however, it is difficult to learn the small
sub-networks directly. The hypothesis is that overparameterized neural
networks are more likely to contain sub-network, which are initialized in such
a way that they can be effectively optimized to solve the task. If finding a
well-initialized sub-network corresponds to winning the lottery then bigger
networks simply have more tickets and are thus more likely to win.

The main positive finding of this paper is that we reproduce the results from
(Najarro and Risi, 2020) with different tasks and a different implementation.
We thus add to the growing body of evidence that shows it is indeed possible
to solve complex tasks with randomly initialized neural networks that learn
using local Hebbian learning rules.

One interesting idea for future work is to re-define an individual
$z=\\{h,k\\}$, i.e. the sampled learning rules $h$ _and_ the learning rule
clusters they are sampled from $k$, such that
$p(z|\theta)=p(h,k|\theta)=\prod_{i=1}^{N}\mathcal{N}(h_{i}|\mu_{k},\sigma)p(k_{i}|\lambda_{i})$.
Without the inner sum, this setup would be much cheaper to optimize for many
learning rule clusters. It might also result in an easier optimization
landscape. However, it would technically not be correct since an individual’s
fitness only depends on the sampled learning rules $h$, not which learning
rule cluster it is sampled from $k$.

A limitation of our approach to Hebbian meta-learning is that the agent cannot
adapt to changes in the reward function during its lifetime, since it does not
observe the reward. The reward is only observed at the evolutionary time
scale. As such all the training and test tasks must share the same reward
function. This is in contrast to common benchmarks in meta RL where the reward
function differs across training and test MDPs. Extending our method with
approaches that can modulate Hebbian learning based on the reward (Abbott,
1990; Krichmar, 2008; Soltoggio et al., 2007, 2018) or indirectly encode
plasticity (Risi and Stanley, 2010), could be a promising future research
direction.

Changes to the pre-registration proposal. The pre-registration proposal
sections (1, 2, 3, and 4) have been changed insignificantly to correct grammar
and spelling. The main findings have been added to the abstract.

## References

  * Abbott (1990) LF Abbott.  Modulation of function and gated learning in a network memory.  _Proceedings of the National Academy of Sciences_ , 87(23):9241–9245, 1990.
  * Andrychowicz et al. (2016) Marcin Andrychowicz, Misha Denil, Sergio Gomez, Matthew W Hoffman, David Pfau, Tom Schaul, Brendan Shillingford, and Nando De Freitas.  Learning to learn by gradient descent by gradient descent.  In _Advances in neural information processing systems_ , pages 3981–3989, 2016.
  * Ba et al. (2016) Jimmy Ba, Geoffrey Hinton, Volodymyr Mnih, Joel Z. Leibo, and Catalin Ionescu.  Using Fast Weights to Attend to the Recent Past.  _ArXiv e-prints_ , Oct 2016.  URL https://arxiv.org/abs/1610.06258v3.
  * Berner et al. (2019) Christopher Berner, Greg Brockman, Brooke Chan, Vicki Cheung, Przemysław Dębiak, Christy Dennison, David Farhi, Quirin Fischer, Shariq Hashme, Chris Hesse, et al.  Dota 2 with large scale deep reinforcement learning.  _arXiv preprint arXiv:1912.06680_ , 2019.
  * Caporale and Dan (2008) Natalia Caporale and Yang Dan.  Spike timing–dependent plasticity: a hebbian learning rule.  _Annu. Rev. Neurosci._ , 31:25–46, 2008.
  * Dempster et al. (1977) Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Duan et al. (2016) Yan Duan, John Schulman, Xi Chen, Peter L Bartlett, Ilya Sutskever, and Pieter Abbeel.  Rl2: Fast reinforcement learning via slow reinforcement learning.  _arXiv preprint arXiv:1611.02779_ , 2016.
  * Finn et al. (2017) Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  _arXiv preprint arXiv:1703.03400_ , 2017.
  * Frankle and Carbin (2018) Jonathan Frankle and Michael Carbin.  The lottery ticket hypothesis: Finding sparse, trainable neural networks.  _arXiv preprint arXiv:1803.03635_ , 2018.
  * Hebb (2005) Donald Olding Hebb.  _The organization of behavior: A neuropsychological theory_.  Psychology Press, 2005.
  * Hinton and Plaut (1987) Geoffrey E Hinton and David C Plaut.  Using fast weights to deblur old memories.  In _Proceedings of the ninth annual conference of the Cognitive Science Society_ , pages 177–186, 1987.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Justesen et al. (2018) Niels Justesen, Ruben Rodriguez Torrado, Philip Bontrager, Ahmed Khalifa, Julian Togelius, and Sebastian Risi.  Illuminating generalization in deep reinforcement learning through procedural level generation.  _arXiv preprint arXiv:1806.10729_ , 2018.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Koch et al. (2015) Gregory Koch, Richard Zemel, and Ruslan Salakhutdinov.  Siamese neural networks for one-shot image recognition.  In _ICML deep learning workshop_ , volume 2. Lille, 2015.
  * Krichmar (2008) Jeffrey L Krichmar.  The neuromodulatory system: a framework for survival and adaptive behavior in a challenging world.  _Adaptive Behavior_ , 16(6):385–399, 2008.
  * Li et al. (2017) Zhenguo Li, Fengwei Zhou, Fei Chen, and Hang Li.  Meta-sgd: Learning to learn quickly for few-shot learning.  _arXiv preprint arXiv:1707.09835_ , 2017.
  * Miconi et al. (2018) Thomas Miconi, Jeff Clune, and Kenneth O Stanley.  Differentiable plasticity: training plastic neural networks with backpropagation.  _arXiv preprint arXiv:1804.02464_ , 2018.
  * Mishra et al. (2017) Nikhil Mishra, Mostafa Rohaninejad, Xi Chen, and Pieter Abbeel.  A simple neural attentive meta-learner.  _arXiv preprint arXiv:1707.03141_ , 2017.
  * Munkhdalai and Trischler (2018) Tsendsuren Munkhdalai and Adam Trischler.  Metalearning with Hebbian Fast Weights.  _ArXiv e-prints_ , Jul 2018.  URL https://arxiv.org/abs/1807.05076.
  * Munkhdalai and Yu (2017) Tsendsuren Munkhdalai and Hong Yu.  Meta Networks.  _ArXiv e-prints_ , Mar 2017.  URL https://arxiv.org/abs/1703.00837v2.
  * Najarro and Risi (2020) Elias Najarro and Sebastian Risi.  Meta-learning through hebbian plasticity in random networks.  _Neural information processing systems (NeurIPS)_ , 2020.
  * Nichol et al. (2018) Alex Nichol, Joshua Achiam, and John Schulman.  On first-order meta-learning algorithms.  _arXiv preprint arXiv:1803.02999_ , 2018.
  * Palm (2020) Rasmus Berg Palm.  Evostrat.  https://github.com/rasmusbergpalm/evostrat, 2020.
  * Pedersen and Risi (2021) Joachim Winther Pedersen and Sebastian Risi.  Evolving and merging hebbian learning rules: Increasing generalization by decreasing the number of rules.  _arXiv preprint arXiv:2104.07959_ , 2021.
  * Prezioso et al. (2018) M Prezioso, MR Mahmoodi, F Merrikh Bayat, H Nili, H Kim, A Vincent, and DB Strukov.  Spike-timing-dependent plasticity learning of coincidence detection with passively integrated memristive circuits.  _Nature communications_ , 9(1):1–8, 2018.
  * Risi and Stanley (2010) Sebastian Risi and Kenneth O Stanley.  Indirectly encoding neural plasticity as a pattern of local rules.  In _International Conference on Simulation of Adaptive Behavior_ , pages 533–543. Springer, 2010.
  * Salimans et al. (2017) Tim Salimans, Jonathan Ho, Xi Chen, Szymon Sidor, and Ilya Sutskever.  Evolution strategies as a scalable alternative to reinforcement learning.  _arXiv preprint arXiv:1703.03864_ , 2017.
  * Santoro et al. (2016) Adam Santoro, Sergey Bartunov, Matthew Botvinick, Daan Wierstra, and Timothy Lillicrap.  Meta-learning with memory-augmented neural networks.  In _International conference on machine learning_ , pages 1842–1850, 2016.
  * Schmidhuber (1992) J. Schmidhuber.  Learning to control fast-weight memories: An alternative to dynamic recurrent networks.  _Neural Computation_ , 4(1):131–139, 1992.
  * Silver et al. (2016) David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al.  Mastering the game of go with deep neural networks and tree search.  _nature_ , 529(7587):484–489, 2016.
  * Snell et al. (2017) Jake Snell, Kevin Swersky, and Richard Zemel.  Prototypical networks for few-shot learning.  In _Advances in neural information processing systems_ , pages 4077–4087, 2017.
  * Soltoggio et al. (2007) Andrea Soltoggio, Peter Durr, Claudio Mattiussi, and Dario Floreano.  Evolving neuromodulatory topologies for reinforcement learning-like problems.  In _2007 IEEE Congress on Evolutionary Computation_ , pages 2471–2478. IEEE, 2007.
  * Soltoggio et al. (2018) Andrea Soltoggio, Kenneth O Stanley, and Sebastian Risi.  Born to learn: the inspiration, progress, and future of evolved plastic artificial neural networks.  _Neural Networks_ , 108:48–67, 2018.
  * Song et al. (2019) Xingyou Song, Wenbo Gao, Yuxiang Yang, Krzysztof Choromanski, Aldo Pacchiano, and Yunhao Tang.  Es-maml: Simple hessian-free meta learning.  _arXiv preprint arXiv:1910.01215_ , 2019.
  * Vinyals et al. (2019) Oriol Vinyals, Igor Babuschkin, Wojciech M Czarnecki, Michaël Mathieu, Andrew Dudzik, Junyoung Chung, David H Choi, Richard Powell, Timo Ewalds, Petko Georgiev, et al.  Grandmaster level in starcraft ii using multi-agent reinforcement learning.  _Nature_ , 575(7782):350–354, 2019.
  * Wang et al. (2016) Jane X Wang, Zeb Kurth-Nelson, Dhruva Tirumala, Hubert Soyer, Joel Z Leibo, Remi Munos, Charles Blundell, Dharshan Kumaran, and Matt Botvinick.  Learning to reinforcement learn.  _arXiv preprint arXiv:1611.05763_ , 2016.
  * Weng (2018) Lilian Weng.  Meta-learning: Learning to learn fast.  _lilianweng.github.io/lil-log_ , 2018.  URL http://lilianweng.github.io/lil-log/2018/11/29/meta-learning.html.
  * Zador (2019) Anthony M Zador.  A critique of pure learning and what artificial neural networks can learn from animal brains.  _Nature communications_ , 10(1):1–7, 2019.