# Second harmonic generation as a minimal model of turbulence

N. Vladimirova1,2 M. Shavit2 S. Belan3,4 G. Falkovich2 1 Brown University,
Providence, RI 02912, USA  2 Weizmann Institute of Science, Rehovot 76100
Israel  3 Landau Institute for Theoretical Physics, 142432 Chernogolovka,
Russia  4 National Research University Higher School of Economics, 101000
Moscow, Russia

###### Abstract

When two resonantly interacting modes are in contact with a thermostat, their
statistics is exactly Gaussian and the modes are statistically independent
despite strong interaction. Considering noise-driven system, we show that when
one mode is pumped and another dissipates, the statistics (of such cascades)
is never close to Gaussian no matter the interaction/noise relation. One finds
substantial phase correlation in the limit of strong interaction (weak noise).
Surprisingly, for both cascades, the mutual information between modes
increases and entropy further decreases when interaction strength decreases.
We use the model to elucidate the fundamental problem of far-from equilibrium
physics: where the information (entropy deficit) is encoded and how singular
measures form. For an instability-driven system (a laser), even a small added
noise leads to large fluctuations of the relative phase near the stability
threshold, while far from it we show that the conversion into the second
harmonic is weakly affected by noise.

##  I Introduction

Second harmonic generation is the simplest fundamental process of nonlinear
wave physics, which is also in the center of numerous practical applications
in laser physics and beyond. The dynamics of the process has been studied
exhaustively SG1 ; SG2 , which cannot be said about statistics, despite the
fact that understanding the influence of noise on the energy conversion is of
paramount practical importance, recently enhanced by the use of meta-materials
SGH1 . Here we address this problem by studying theoretically a two-mode
resonant system driven by a combination of pumping and random noise. Our
motivation is two-fold. Apart from the classical conversion problem, we find
this system ideally suited for elucidating the fundamental problems of non-
equilibrium physics. When one mode is stochastically forced and another is
dissipated, that presents a minimal model of turbulence cascade. The freedom
to force either mode allows us to elucidate the basic differences between
direct and inverse cascades. Apart from energy, we shall be interested in the
entropy of such far-from equilibrium state, which is expected to be much lower
than in thermal equilibrium with the same energy.

A remarkable property, common for all systems of resonantly interacting waves
and shared with hydrodynamic systems VSF , is that the canonical thermal
equilibrium has exactly Gaussian statistics, and the modes fluctuate
independently, regardless of the interaction strength. Here we describe how
deviations from equilibrium diminish entropy and build correlations between
the two modes. Far from equilibrium the joint two-mode statistics is never
close to Gaussian, even when the marginal distribution of every mode is close
to Gaussian. On the one hand, the entropy decrease means that the statistical
distribution is getting more non-uniform, which poses the question: Can it
lead all the way to singularity like the celebrated Sinai-Ruelle-Bowen
measures in dynamical systems D ; FS ? We show that this is indeed so: the
measure in the phase space is getting singular in the double limit of strong
non-equilibrium and weak interaction. On the other hand, since entropy is
missing information, any entropy decrease poses another question: Where all
this extra information about non-equilibrium is encoded? First, we find out
how the entropies of the three marginal distributions, of each mode amplitude
and their phase difference, go down as the system deviates from equilibrium.
Second, we find out which part of the entropy decrease is due to inter-mode
correlation. This is properly measured by the mutual information (rather than
by the pair correlation function, suitable for Gaussian statistics only).

The process of the second harmonic generation is described by the following
model Hamiltonian (assuming perfect resonance)

$\displaystyle{\cal
H}_{0}=\omega|a_{1}|^{2}+2\omega|a_{2}|^{2}+Va_{1}^{*2}a_{2}+V^{*}a_{1}^{2}a_{2}^{*}\
.$ (1)

Here $a_{1}$ and $a_{2}$ are the complex amplitudes of two non-linearly
coupled modes having frequencies $\omega$ and $2\omega$, respectively, and $V$
is the interaction constant (considered real positive without loss of
generality). The two coupled complex equations govern dynamics:
$\dot{a}_{k}=-i{\partial{\cal H}_{0}}/{\partial a_{k}^{*}}$, $k=1,2$. We
eliminate the linear terms in these equations by introducing the envelopes

$b_{1}=a_{1}e^{i\omega t},b_{2}=a_{2}e^{2i\omega t}.$ (2)

That results in a strongly interacting system with a cubic Hamiltonian ${\cal
H}=Vb_{1}^{*2}b_{2}+V^{*}b_{1}^{2}b_{2}^{*}$.

Due to the symmetry $b_{1}\to b_{1}e^{i\phi}$, $b_{2}\to b_{2}e^{2i\phi}$ the
system $\dot{b}_{k}=-i{\partial{\cal H}}/{\partial b_{k}^{*}}$ has an extra
integral of motion $N=|b_{1}|^{2}+2|b_{2}|^{2}$ and is completely integrable;
the phase portrait is presented in Appendix A.1. Let us add dissipation and
stochastic pumping:

$\displaystyle\dot{b}_{1}=-2iV^{*}b_{1}^{*}b_{2}-\gamma_{1}b_{1}+\xi_{1}(t),$
(3) $\displaystyle\dot{b}_{2}=-iVb_{1}^{2}-\gamma_{2}b_{2}+\xi_{2}(t).$ (4)

Here $\gamma_{1}$ and $\gamma_{2}$ are the damping coefficients, and $\xi_{1}$
and $\xi_{2}$ are independent Gaussian random forces with zero mean
$\langle\xi_{i}(t)\rangle=0$ and the variance
$\langle\xi_{i}(t_{1})\xi_{j}^{*}(t_{2})\rangle=P_{i}\delta_{ij}\delta(t_{2}-t_{1})$.

We mainly focus on the properties of the statistically steady solutions of the
system (3,4) in the case when one mode is forced, while the other is damped.
Since the modes enter the Hamiltonian in a non-symmetric way, there are two
possibilities: one either can pump the first (lower frequency) mode and dump
the second (higher frequency) mode or vice versa. The former scenarios
qualitatively corresponds to the direct energy cascade, while the second is
reminiscent to the inverse cascade.

We wish to understand how much information is needed in order to build a
turbulent state and how much one learns about one mode by observing another.
For that we will use the metrics from information theory: entropy and mutual
information. The answer to the first question is given by the decrease in
entropies

$\displaystyle S_{12}=-\int
db_{1}db_{1}^{*}db_{2}db_{2}^{*}\,\rho(b_{1},b_{2})\ln\rho(b_{1},b_{2})\,,$
$\displaystyle S_{1}=-\int db_{1}db_{1}^{*}\,\rho(b_{1})\ln\rho(b_{1})\,,$
$\displaystyle S_{2}=-\int db_{2}db_{2}^{*}\,\rho(b_{2})\ln\rho(b_{2})\,,$

where $\rho$ is either full or marginal probability distribution. The answer
to the second question is given by the mutual information between the modes:

$I_{12}=S_{1}+S_{2}-S_{12}\ .$ (5)

Fig.1 demonstrates the growth of the mutual information versus the degree of
non-equilibrium (an analog of the Reynolds number defined below, see (8)).

As one of the simplest model of energy transfer, the system of two coupled
oscillators has received considerable attention in the literature Kumar_2008 ;
Bonetto_2004 ; Ciliberto_2013 ; Falasco_2015 ; Chun_2015 ; Mura_2018 ;
Horowitz_2019 ; Courant1 ; Courant . In particular, in the mathematical
literature, one finds an analysis of a two-mode system with a quadratic
Hamiltonian ${\cal H}=Ta_{1}^{*2}a_{2}^{2}$ with the purpose to get insight
into the energy transfer in wave turbulence Courant1 ; Courant . What
distinguishes our model is that it directly corresponds to physical reality
and allows experimental validation. In addition, an asymmetry between the
modes allows us to compare direct and inverse cascades, which turn out quite
different. Another distinction is that we add entropic and informational
consideration to the energetic analysis.

Figure 1:  Mutual information versus Reynolds number for direct and inverse
cascades (red and blue lines respectively). Three dimensional distribution are
computed with bin size $\Delta\theta=2\pi/32$ and
$\Delta\rho^{2}_{1,2}/n_{1,2}=1$ (circles) and 0.5 (diamonds).

The paper is organized from the point of view of entropy: we start from the
maximal-entropy equilibrium and investigate near-equilibrium states in Section
II. We then move to study the noise-driven direct and inverse cascades in
Section III. We define a dimensionless parameter akin to the Reynolds number
and show how entropy decreases as the Reynolds number increases. We begin with
the limit of small Reynolds number in III.1 and continue to large Reynolds
number in III.2, where entropy decreases all the way down as the measure
becomes singular in the limit of infinite Reynolds number. In Section IV we
consider an instability-driven first harmonic and study the influence of noise
on the conversion process, this can serve as a simple model for a laser
generating second harmonic. Conclusion V briefly lists our main results.

##  II Near Thermal equilibrium

Adopting the language of stochastic thermodynamics, one can call the ratios
$P_{1}/\gamma_{1}\equiv T_{1}$ and $2P_{2}/\gamma_{2}\equiv T_{2}$ effective
temperatures experienced by two modes $b_{1}$ and $b_{2}$ which are governed
by the Langevin equations (3) and (4). If $\Delta T=T_{1}-T_{2}=0$, then it is
straightforward to find from the Fokker-Planck equation or entropy maximum the
steady-state joint probability distribution:

${\cal P}_{0}={1\over
Z}\exp\biggl{(}-\frac{2|b_{1}|^{2}+4|b_{2}|^{2}}{T}\biggr{)}={1\over
Z}\exp\biggl{(}-\frac{2N}{T}\biggr{)}.$ (6)

Despite strong interaction, this distribution is exactly Gaussian and the
modes are statistically independent. The later means that the mean energy flux
and the mutual information between modes are zero. Thermal equilibrium
corresponds to the equipartition of the quadratic invariant:
$\langle|b_{1}|^{2}\rangle\equiv n_{1}=2n_{2}\equiv
2\langle|b_{2}|^{2}\rangle$.

What can we say about the system’s statistics when modes are subject to
different effective temperatures? Let us introduce the dimensionless measure
of non-equilibrium

$\sigma=\frac{\Delta T}{T},\ \ \text{where}\ \
T=\frac{P_{1}+2P_{2}}{2(\gamma_{1}+\gamma_{2})}.$ (7)

Another dimensionless parameter quantifies interaction strength relative to
the dissipation:

$\chi={(\gamma_{1}+\gamma_{2})^{3}\over(P_{1}+2P_{2})|V|^{2}}\ .$ (8)

Figure 2:  Top row: deviation of entropies from equilibrium ($\Delta T=0$) and
mutual information for $\gamma_{1}=\gamma_{2}$ and different strength of
interaction marked by color. Equilibrium is Gaussian with $S_{{\rm
eq}1}=1/\ln(2)$, $S_{{\rm eq}2}=1/\ln(2)-1$, and $S_{{\rm
eq}\theta}=\log_{2}(2\pi)$. Dashed lines show the entropy deviation of a
single mode respectively at $T_{1}$ and $T_{2}$ for marginal distributions,
$\Delta S_{1}=\log_{2}(1+\Delta T/4T)$ and $\Delta S_{2}=\log_{2}(1-\Delta
T/4T)$ . The entropies of individual modes are affected by the change of
variance of (nearly Gaussian) distribution and by deviation from the Gaussian
shape. For close-to-equilibrium cases presented here the first effect is
dominant. Bottom row: panels (e) and (f) show deviation of entropies from the
entropies of Gaussian distributions with the same variance; panel (g)
illustrates the quadratic dependence of the entropy of phase distribution,
panel (h) illustrates the quadratic dependence for mutual information (solid
lines for $\Delta T>0$; dashed lines for $\Delta T<0$). The offset 0.005 is
the effect of finite ensemble and bin size; the offset decreases as the size
of data set and resolution improve. 3D distributions are computed with bin
sizes $2\pi/32$ for the phase and $0.1T$ for $\rho_{1,2}$.

Denote $\rho_{1,2}=|b_{1,2}|$ and $\theta=\arg(b_{1}^{2}b_{2}^{*})$. From (3)
and (4), the steady-state equations on the second moments read:

$\displaystyle-4V\langle\rho_{1}^{2}\rho_{2}\sin\theta\rangle-2\gamma_{1}\langle\rho_{1}^{2}\rangle+P_{1}=0,$
(9) $\displaystyle
2V\langle\rho_{1}^{2}\rho_{2}\sin\theta\rangle-2\gamma_{2}\langle\rho_{2}^{2}\rangle+P_{2}=0,$
(10)

The time derivative of the real part of the third moment is given by:
$d\langle{\cal H}\rangle/dt=-(2\gamma_{1}+\gamma_{2})\langle\cal H\rangle$,
since
$\langle\xi_{1}b_{1}b_{2}^{*}\rangle=\langle\xi_{2}^{*}b_{1}^{2}\rangle=0$.
Therefore, in any steady state, either in thermal equilibrium or out of it,
one has

$\langle{\cal H}\rangle=2V\langle\rho_{1}^{2}\rho_{2}\cos\theta\rangle=0.$
(11)

Equations (9)-(11) are valid for any values of $\sigma$ and $\chi$.

At $\sigma\neq 0$, the probability density ${\cal P}(b_{1},b_{2})$ is non-
Gaussian in non-equilibrium, yet it is close to Gaussian when $|\sigma|\ll 1$
for all values of $\chi$. The simplest to treat is the limit of small
interaction, $\chi\gg 1$. In this case, the first correction to Eq. (6) is
determined by the energy flux between modes, which is small and proportional
to the temperature difference:

$\displaystyle\ln{\cal P}(b_{1},b_{2})\approx-{2|b_{1}|^{2}\over
T_{1}}-{4|b_{2}|^{2}\over T_{2}}-$ (12) $\displaystyle-{4\Delta T\over
P_{1}T_{2}+P_{2}T_{1}}{\rm\,Im\,}[V^{*}b_{1}^{*2}b_{2}]+O(\chi^{-2})$ (13)

Smallness of interaction multiplies the parameter of non-equilibrium $\Delta
T/T$ in the right hand side of Eq. (12), so that this result is valid even
when $\Delta T/T$ is not small. That means that, as long as both temperatures
remain finite and interaction is weak, even far from equilibrium the relative
entropy is small:

$D({\cal P}|{\cal P}_{0})=\int db_{1}db_{1}^{*}db_{2}db_{2}^{*}{\cal
P}\ln({\cal P}/{\cal P}_{0})\propto\chi^{-1}\ll 1\ ,$

as well as the mutual information.

In the opposite limit $\chi\ll 1$ or $V\to\infty$, the non-Gaussian correction
is again proportional to the product of the degree of non-equilibrium and the
small parameter $\chi$. In terms of $x=|b_{1}|^{2}$ and $y=2|b_{2}|^{2}$ we
obtain:

$\ln{\cal P}\approx-{x+y\over T}+{\chi^{1/2}\Delta T\over T}f(x,y,\theta)\ ,$

where the correction satisfies the equation

${2x{y}\over\sqrt{yT}}\left[\sin\theta\left({\partial\over\partial
y}-{\partial\over\partial x}\right)+{x-2y\over
2xy}\cos\theta{\partial\over\partial\theta}\right]f=x-y\,.$

The correction is odd in phase difference, $f(-\theta)=-f(\theta)$, and scales
linearly with amplitudes, so that it is substantial at small amplitudes. In
the limits, $f\to g(\theta)\sqrt{y/2T}$, where $g=\sin\theta$ at $y\ll x$, and
$g=\int d\theta/\cos\theta$ at $y\gg x$.

It makes sense to compare entropies at the same mean quadratic energy $N$. To
see how entropy goes down on the way to turbulence we shall subtract the total
entropy from its maximal equilibrium value, which quantifies the amount of
information one needs to create a turbulent state: $\Delta
S(N,n_{2}/n_{1})=S_{0}-S_{12}$.

Numerics support quadratic decrease of $S_{12}(\Delta/T)$ and increase of
$I_{12}(\Delta/T)$ up to $\Delta T\simeq 4T$, see Figure 2.

When $\Delta T/T$ exceeds one, the functions are not even which demonstrates
the statistical difference between upward and downward energy conversion. We
see stronger deviations from Gaussianity for negative $\Delta T<0$, which
corresponds to the downward energy flow and to an inverse cascade at $\Delta
T/T\to-\infty$. The physical difference is that the first mode pumps the
second one as an additive force, while the second mode pumps the first one as
a multiplicative instability. Therefore, it seems natural that the entropy is
generally lower and the mutual information higher for an inverse transfer. The
analysis of the separate distributions of two amplitudes and the relative
phase shows that the entropy of the driven mode (say, $S_{1}(\Delta/T)$ for a
direct transfer) grows with $\Delta T$ slower than the entropy of the
dissipated mode and $S_{\theta}(\Delta/T)$ decrease, see panels (g) and (h) in
Figure 2.

##  III Turbulent cascades

Now let us have an energy cascade in our model: pumping one mode and
dissipating another. When energy flows from lower frequency mode to higher,
i.e $0=\xi_{2}=\gamma_{1}$ in Eqs. (3) and (4), the cascade is called direct,
and inverse when $0=\xi_{1}=\gamma_{2}$. In these cases, the only
dimensionless parameter is $\chi=\gamma^{3}/(P|V|^{2})$, where $P$ is the
intensity of noise acting upon the driving mode, and $\gamma$ denotes the
damping coefficient of the dissipating mode. As we shall see below, $\chi$ to
some extent plays the role of the Reynolds number of hydrodynamics in a sense
that it determines how low is the entropy and how much the occupation numbers
deviate from the equipartition $n_{1}=2n_{2}$, even though the system is not
close to thermal equilibrium for however small or large $\chi$.

Balance of the quadratic invariant, $N$, means that the dissipating mode keeps
the magnitude of order of its equilibrium value: $n_{2}=P/4\gamma$ for the
direct cascade and $n_{1}=P/\gamma$ for the inverse cascade. How much the mode
which is pumped exceeds the equipartition value is determined by the value of
$\chi$, as described below. Note that this parameter can be interpreted as the
squared ratio of the dissipation rate $\gamma$ and the nonlinear transfer rate
$Vn\simeq V\sqrt{P/\gamma}$.

When $\chi$ is small, the interaction between modes is strong and the energy
transfer is fast, so that the occupation numbers are expected to be close to
equipartition, yet the statistics is not expected to be close to separable
Gaussian form given by (6). Even though the noise is weak, it is white, that
is a singular perturbation destroying integrability everywhere in the phase
space Kurchan ; we shall see below how non-trivial the probability
distribution is already in this limit.

One may naively expect that in the opposite limit of large $\chi$, when the
noise is strong and interaction is weak, the correlation between modes would
be weak too. We shall show below that the opposite is true far from
equilibrium: the necessity to carry the flux makes the modes strongly
correlated precisely because of a strong noise and weak interaction. It is in
this limit we find the lowest entropy and the maximal mutual information
between modes, as well as appearance of singular measure in phase space.

.

Figure 3:  Probabilities of the occupation numbers and the phase for: (a, b)
inverse cascade, small $\chi$; (c, d) direct cascade, small $\chi$; (e, f)
inverse cascade, large $\chi$; (g, h) direct cascade, large $\chi$. Each
dataset contains 20M datapoints, at $\Delta t=0.01$ for inverse cascade and
large $\chi$ and $\Delta t=0.1$ for all other cases. For the inverse cascade,
$\chi=\gamma_{1}^{3}/(2P_{2}V^{2})$ and $\nu_{1}=P_{2}/\gamma_{1}$. We use
$\gamma_{1,2}=0.01$ for small $\chi$ and $\gamma_{1,2}=1$ for large $\chi$. In
all cases, $V=1$. Broken lines in (g) correspond to the approximation (31,36).

###  III.1 Small Reynolds number: strong-interaction-weak-noise limit

####  III.1.1 Inverse cascade

In an inverse cascade, energy goes from high to low frequency, so we set
$\gamma_{2}=0$ and $P_{1}=0$ in (3,4):

$\displaystyle\dot{b}_{1}=-2iV^{*}b_{1}^{*}b_{2}-\gamma b_{1},$ (14)
$\displaystyle\dot{b}_{2}=-iVb_{1}^{2}+\xi(t).$ (15)

In the steady state, the energy input rate $P$ must be equal to the
dissipation rate $\gamma n_{1}$ and to the energy flux from the second mode to
the first given by the imaginary part of the third cumulant:
$2V\langle\rho_{1}^{2}\rho_{2}\sin\theta\rangle=-P$. Then, from the energy
balance we obtain $n_{1}=P/\gamma$ and
${2V}\langle\rho_{1}^{2}\rho_{2}\sin\theta\rangle=-{P}$, so that
$\langle\rho_{1}^{2}\rho_{2}\rangle\geq{P}/{2V}$. Also, from the condition
$\frac{d}{dt}\langle\ln|b_{1}|^{2}\rangle=0$ we find
$2V\langle{\rho_{2}}\sin\theta\rangle=-{\gamma}$ and, therefore
$\langle{\rho_{2}}\rangle\geq{\gamma}/{2V}$.

From Eqs. (14,15) it is straightforward to see that when $\chi\ll 1$ the
steady-state probability distribution ${\cal P}(b_{1},b_{2})$ cannot be close
to the equilibrium Gaussian (6) with the temperature $T=2P/\gamma$ and the
equipartition $P/\gamma=n_{1}=2n_{2}$. Indeed, the stationarity of
$\langle{\cal
H}^{2}\rangle=4|V|^{2}\langle\rho_{1}^{4}\rho_{2}^{2}\cos^{2}\theta\rangle$
requires $\langle{\cal
H}^{2}\rangle/\langle\rho_{1}^{4}\rangle=\frac{|V|^{2}P}{2\gamma}$, but this
contradicts the Gaussian ratio which gives $\langle{\cal
H}^{2}\rangle/\langle\rho_{1}^{4}\rangle=\frac{\left|V\right|^{2}P}{\gamma}$.

Thus, small value of $\chi$ does not mean that the system is near
equipartition. In agreement with this conclusion, results of numerical
modeling presented at two lower left panels of Figure 3 show that at neither
of marginal distributions of the mode amplitudes is close to Gaussian and that
the phase distribution does not become uniform as $\chi$ tends to zero. This
is also reflected in nonzero value $I_{12}(+0)$ of mutual information between
modes in this limit, see the right panel of Figure 1.

####  III.1.2 Direct cascade

Direct cascade corresponds to the choice $\gamma_{1}=0$ and $P_{2}=0$ in
(3,4):

$\displaystyle\dot{b}_{1}=-2iV^{*}b_{1}^{*}b_{2}+\xi(t),$ (16)
$\displaystyle\dot{b}_{2}=-iVb_{1}^{2}-\gamma b_{2}.$ (17)

Now energy goes from low to high frequency. Again, in the steady-state regime,
the energy input rate $P/2$ must be equal to the dissipation rate $2\gamma
n_{2}$ and to the energy flux from the first mode to the second given by the
imaginary part of the third cumulant:
$2V\langle\rho_{1}^{2}\rho_{2}\sin\theta\rangle=P/2$, so $n_{2}=P/4\gamma$.
From $\frac{d}{dt}\langle\ln|b_{2}|^{2}\rangle=0$ we find
$\langle\rho_{1}^{2}\sin\theta/\rho_{2}\rangle={\gamma}/{V}$. Therefore
$\langle\rho_{1}^{2}\rho_{2}\rangle\geq{P}/{4V}$ and
$\langle{\rho_{1}^{2}}/{\rho_{2}}\rangle\geq{\gamma}/{V}$.

When $\chi\to 0$, the dimensionless flux
$\langle\rho_{1}^{2}\rho_{2}\sin\theta\rangle/n_{1}n_{2}^{1/2}=\chi^{1/2}$ is
small, which may suggest that phase-space distribution is close to the
Gaussian equilibrium (6) with $T=P/\gamma$ and that the phase distribution is
close to uniform. Furthermore, as opposed to the case of inverse cascade
discussed above, the equality obtained from the stationarity of $\langle{\cal
H}^{2}\rangle$:

$\langle{\cal
H}^{2}\rangle=4|V|^{2}\langle\rho_{1}^{4}\rho_{2}^{2}\cos^{2}\theta\rangle=\frac{2|V|^{2}P}{\gamma}\langle\rho_{1}^{2}\rho_{2}^{2}\rangle\
,$ (18)

is achieved by the Gaussian distribution with $n_{1}=2n_{2}=P/2\gamma$.
However, numerical data, as can be seen from two upper left panels of Figure
3, shows that even though the marginal distributions of amplitudes are close
to Gaussian with equipartition, $n_{1}\approx 2n_{2}$, the phase distribution
is far from flat and deviation from equilibrium is substantial. The mutual
information between modes as a function of $\chi$ exhibits a non-zero value of
$I_{12}(+0)$ (see Figure 1) which is also a clear footprint of non-
equilibrium.

###  III.2 Large Reynolds number: weak interaction, strong noise limit

####  III.2.1 Inverse cascade

The pair of complex equations (14,15) can be rewritten as three real ones
since the overall phase drops out:

$\displaystyle\dot{\rho}_{1}=-2|V|\rho_{1}\rho_{2}\sin\theta-\gamma\rho_{1},$
(20) $\displaystyle\dot{\rho}_{2}=|V|\rho_{1}^{2}\sin\theta+{P\over
4\rho_{2}}+{\zeta(t)\over\sqrt{2}},$
$\displaystyle\dot{\theta}=|V|\frac{\rho_{1}^{2}-4\rho_{2}^{2}}{\rho_{2}}\cos\theta+{\frac{\zeta(t)}{\sqrt{2}\rho_{2}}},$
(21)

where $\zeta(t)$ is the real white noise with zero mean
$\langle\zeta(t)\rangle=0$ and the pair correlation function
$\langle\zeta(t_{1})\zeta(t_{2})\rangle=P\delta(t_{1}-t_{2})$.

When $\chi\gg 1$, Eqs. (20)-(21), can be further simplified by assuming that
relative phase is locked on $\theta=-\pi/2$ most the time. Then, one gets the
following closed equations for the amplitudes dynamics

$\displaystyle\dot{\rho}_{1}=2|V|\rho_{1}\rho_{2}-\gamma\rho_{1},$ (22)
$\displaystyle\dot{\rho}_{2}=-|V|\rho_{1}^{2}+{P\over
4\rho_{2}}+{\zeta(t)\over\sqrt{2}}.$ (23)

A hypothesis that the modes are statistically independent in this limit is
shown incorrect in the Appendix A.4. This result is in sharp contrast with the
model described in Courant1 , where authors found the factorized joint
probability density ${\cal P}(\rho_{1},\rho_{2})$ of mode amplitudes in the
limit when their analogue of the parameter $\chi$ is large.

While constructing the probability densities for inverse cascade at
$\chi\to\infty$ turns out to be a tricky task, it is straightforward to
describe general features of stochastic dynamics dictated by Eqs. (22) and
(23). Namely, these pair of nonlinearly coupled equations suggest the
following cyclical evolution: $\rho_{1}$ stays close to zero most of the time
while $\rho_{2}$ undergoes diffusion in a repulsive logarithmic potential;
when $\rho_{2}$ sufficiently outgrows the threshold level $\gamma/2|V|$,
$\rho_{1}$ shoots up and quickly diminishes $\rho_{2}$; after that $\rho_{1}$
also resets to the near-zero level and the stochastic dynamics of $\rho_{2}$
starts from scratch. The mode dynamics during the intermittent burst events
can be described by simplified equations

$\displaystyle\dot{\rho}_{1}=2|V|\rho_{1}\rho_{2}-\gamma\rho_{1},$ (24)
$\displaystyle\dot{\rho}_{2}=-|V|\rho_{1}^{2}.$ (25)

Compared with Eqs. (22) and (23), we neglected the terms associated with
noise.

Figure 4: An individual realization of the modes trajectories during one of
the events. The continues lines are obtained from numerical simulations of
Eqs. (14) and (15, and the dashed lines represent theoretical fit based on
Eqs. (26) and (27)).

Equations (24) and (25) are exactly solvable yielding

${\rho}_{1}^{2}(t)={r_{1}^{2}-2\left(\rho_{2}(t)-\frac{\gamma}{2V}\right)^{2}+2\left(r_{2}-\frac{\gamma}{2|V|}\right)^{2}}\,,$
(26)
$\displaystyle{\rho}_{2}(t)=\frac{\gamma}{2V}+\sqrt{\left(r_{2}-\frac{\gamma}{2V}\right)^{2}+\frac{r^{2}_{1}}{2}}\tanh\left[\frac{1}{2}\ln\frac{\sqrt{\bigl{(}r_{2}-\frac{\gamma}{2V}\bigr{)}^{2}+\frac{r^{2}_{1}}{2}}+r_{2}-\frac{\gamma}{2V}}{\sqrt{\bigl{(}r_{2}-\frac{\gamma}{2V}\bigr{)}^{2}+\frac{r^{2}_{1}}{2}}-r_{2}+\frac{\gamma}{2V}}-2Vt\sqrt{\Bigl{(}r_{2}-\frac{\gamma}{2V}\Bigr{)}^{2}+\frac{r^{2}_{1}}{2}}\right],$
(27)

where $r_{1}=\rho_{1}(0)$, $r_{2}=\rho_{2}(0)$ are the initial conditions.
Estimating $r_{2}\sim{\gamma/V}$ and $r_{1}\ll r_{2}$, we see from Eq. (27)
that the duration of such burst event is $\sim\gamma^{-1}$, which is much
smaller than the typical inter-events period $\sim{\gamma^{2}}/{PV^{2}}$. As
follows from Eq. (26), the amplitude of the dissipated mode grows from the
initial value $r_{1}\ll r_{2}$ to the maximum value
$\rho_{1\text{max}}=\sqrt{r_{1}^{2}+2(r_{2}-\frac{\gamma}{2V})^{2}}\approx\sqrt{2}(r_{2}-\frac{\gamma}{2V})$
(attaining it at the moment when $\rho_{2}(t)=\frac{\gamma}{2V}$) and finally
returns to the starting level $\rho_{1}(\infty)=r_{1}$. Such bursts are likely
responsible for pulses running in shell models, which are chains of
interacting triplets.

Numerical simulations confirm the intermittent nature of system dynamics
described above. Namelly, panel (f) of Figure 3 reveals that the relative
phase is indeed locked at $-\pi/2$. Figure 4 illustrates that analytical
prediction based on the assumption of phase locking (see Eqs. (26) and (27))
are in excellent agreement with numerical data extracted from simulations of
Eqs. (14) and (15). As visible in panel (e) of Figure 3, the tails of the
amplitudes probability distributions strongly depend on $\chi$; the fits of
${\cal P}(\rho_{1})$ and ${\cal P}(\rho_{2})$ by the broken lines in Figure 3
are empirical. As for the mutual information, from Figure 1 we see that
$I_{12}\propto\ln\chi$ for inverse cascade with $\chi\gg 1$.

####  III.2.2 Direct cascade

We conclude the treatment of our turbulent cascades with the case of a direct
cascade in the limit of large Reynolds number, $\chi\rightarrow\infty$. It is
in this limit we find a window to the way a singular measure is formed far
away from equilibrium. We indeed find that in this limit the full probability
distribution is singular and thus corresponds to the lowest entropy state
$S_{12}\rightarrow-\infty$.

In addressing the weak-interaction-strong-noise limit $\gamma\to\infty$
($\chi\rightarrow\infty$), it is convenient to express $b_{2}$ from (17) as an
integral, which in the leading order shows that the second amplitude is
enslaved to the first one:

$b_{2}(t)=-iV\int_{-\infty}^{t}b_{1}^{2}(t^{\prime})e^{\gamma(t^{\prime}-t)}dt^{\prime}\approx-\frac{iVb_{1}^{2}}{\gamma}.$
(28)

Substituting this relation into (16), we get a closed equation on the
stochastic dynamics of $b_{1}$

$\dot{b}_{1}=-\frac{2V^{2}}{\gamma}b_{1}^{*}b_{1}^{2}+\xi(t),$ (29)

from which one finds the following expressions for the marginal probability
distributions

$\displaystyle{\cal P}(b_{1})$
$\displaystyle=Z_{1}^{-1}\exp\biggl{(}-\frac{2V^{2}}{\gamma
P}|b_{1}|^{4}\biggr{)},$ (30) $\displaystyle{\cal P}(b_{2})$
$\displaystyle=Z_{2}^{-1}\exp\biggl{(}-\frac{2\gamma}{P}|b_{2}|^{2}\biggr{)},$
(31)

which are valid at $\rho_{1},\rho_{2}\ll{\gamma}/{|V|}$. Thus, the whole
probability density in the four-dimensional phase space is singular at
$\chi\to\infty$, sitting on a three-dimensional manifold

$\\!\\!\\!{\cal P}(b_{1},b_{2})=\frac{2}{\sqrt{\pi
P\gamma}}\exp\biggl{(}-\frac{4V^{2}|b_{1}|^{4}}{P\gamma}\biggr{)}\delta(b_{2}-\frac{iV}{\gamma}b_{1}^{2})\,,$
(32)

so that the total entropy $S_{12}\to-\infty$. Note also that Eq. (32) yields
large ratio of the typical mode amplitudes:
$\rho_{1}^{2}/\rho_{2}^{2}\simeq\sqrt{\gamma^{3}/P|V|^{2}}=\sqrt{\chi}\gg 1$.

Since the distribution over the overall phase is flat, one can integrate it
out and conclude that the distribution in the the three-dimensional space of
variables $\rho_{1},\rho_{2},\theta$ concentrates on the curve
$\rho_{2}\propto\rho_{1}^{2}$. Interestingly, with increasing $\chi$ the joint
probability distribution ${\cal P}(b_{1},b_{2})$ is getting sharper than
Gaussian along this curve. This is different from the model discussed in
Courant1  where the driving mode is nearly Gaussian in this limit, the
relative amplitudes of both modes fluctuate, and only the relative phase is
fixed, that is the joint probability density is only singular with respect to
the phase difference $\theta$.

Considering large but finite $\chi$, instead of delta-function in Eq. (32) one
obtains the distribution with a finite width which is the variance of the
difference $b_{2}-Vb_{1}^{2}/\gamma$. To estimate this width we further expand
Eq. (28)

$\displaystyle b_{2}(t)\approx$ $\displaystyle-
iV\int_{-\infty}^{t}\left[b_{1}^{2}(t)+(t^{\prime}-t){db_{1}^{2}(t)\over
dt}\right]e^{\gamma(t^{\prime}-t)}dt^{\prime}$ (33) $\displaystyle=$
$\displaystyle-{iVb_{1}^{2}(t)\over\gamma}+{2iVb_{1}(t)\over\gamma^{2}}{db_{1}(t)\over
dt}.$ (34)

From Eq. (33) we get

$\langle|b_{2}+iVb_{1}^{2}/\gamma|^{2}\rangle={4|V|^{2}\over\gamma^{4}}\langle|b_{1}\dot{b}_{1}|^{2}\rangle={4|V|^{2}P^{2}\over\gamma^{4}}\
.$ (35)

Dividing this result by $\langle|b_{2}|^{2}\rangle=P/\gamma$, one obtains that
the relative squared width behaves as $1/\chi$ (as expected, it tends to zero
when $\chi\to\infty$). The non-zero width at finite values of $\chi$ entails
the finite entropy of the distribution ${\cal P}(b_{1},b_{2})$:
$S_{12}\simeq-\ln\chi$. This analytical prediction is supported by Figure 1.
For comparison, the naive Gaussian ansatz yields $S_{1}+S_{2}\simeq\ln
n_{1}n_{2}\propto\ln\chi^{-1/2}$, since $n_{2}=P/4\gamma$ and
$n_{1}\simeq\sqrt{P\gamma}/|V|$.

Now let us plug Eq. (28) into Eq.(17), then solving the corresponding steady
Fokker-Planck equation (see A.3) one arrives at the next order correction for
the marginal probability distribution for $\rho_{1}\ll{\gamma}/{|V|}$:

${\cal
P}(\rho_{1})\propto\rho_{1}(1-\frac{4|V|^{2}}{\gamma^{2}}\rho_{1}^{2})\exp(-\frac{2|V|^{2}}{\gamma
P}\rho_{1}^{4}+\frac{16|V|^{4}}{3\gamma^{3}P}\rho_{1}^{6})\,.$ (36)

which is more accurate than Eq. (30). Unfortunately, extracting similar
correction for the probability density ${\cal P}(\rho_{2})$ as well as the
further high order corrections to ${\cal P}(\rho_{1})$ is more challenging.

As can be seen from panel (g) of Figure 3, Eqs. (31) and (36) allow us to fit
the numerical data quite accurately in the range
$\rho_{1},\rho_{2}\ll{\gamma}/{|V|}$. Expectedly, the agreement between
numeric and analytic results improves with the growth of $\chi$.

##  IV Laser generation

Let us now pump the first harmonic by an instability, for instance, in a
laser. Consider first pure dynamics, setting $\xi_{1}=\xi_{2}=0$ in (3,4) and
changing sign in front of $\gamma_{1}$, which now describes gain for an
optical signal. Then the resulting evolution satisfies three closed equations:

$\displaystyle{d\rho_{1}^{2}\over dt}=-2F+2\gamma_{1}\rho_{1}^{2},$ (37)
$\displaystyle{d\rho_{2}^{2}\over dt}=F-2\gamma_{2}\rho_{2}^{2},$ (38)
$\displaystyle{dF\over
dt}=(2\gamma_{1}-\gamma_{2})F+2|V|^{2}\bigl{(}\rho_{1}^{4}-4\rho_{1}^{2}\rho_{2}^{2}\bigr{)},$
(39)

where $F=2V\rho_{1}^{2}\rho_{2}\sin\theta$ is the instantaneous flux.

Apart from the trivial unstable fixed point $\rho_{1}=\rho_{2}=0$, Eqs.
(37)-(39) have the stationary point
$\bar{\rho}_{1}^{2}=\gamma_{1}\gamma_{2}/2|V|^{2}$,
$\bar{\rho}_{2}^{2}=\gamma_{1}^{2}/4|V|^{2}$, $\bar{\theta}=\pi/2$ (and, thus,
$\bar{F}=\gamma_{2}\gamma_{1}^{2}/2|V|^{2}$). This means that in the the
degenerate case $\gamma_{2}=2\gamma_{1}$, the system possesses the steady
state $\rho_{1}=2\rho_{2}$ existing for any $\theta$. This marginal stability
turns into an instability of the steady state at $\gamma_{2}<2\gamma_{1}$ and
into stability at $\gamma_{2}>2\gamma_{1}$. In what follows, we consider
$\gamma_{2}>2\gamma_{1}$. Note that in practice both $\gamma_{1}$ and
$\gamma_{2}$ often depend on the amplitudes, for instance, due to gain
saturation or/and nonlinear damping. However, our main focus here is on the
noise impact on the steady state, so we will treat $\gamma_{1},\gamma_{2}$
taken near this state as constants.

Let us now add a random pumping and study its influence on the efficiency and
statistics of conversion. The modes amplitudes, $\rho_{1}$ and $\rho_{2}$ and
the relative phase $\theta$ are governed by the following equations

$\displaystyle\dot{\rho}_{1}=-2V\rho_{1}\rho_{2}\sin\theta+\gamma_{1}\rho_{1}+{P\over
4\rho_{1}}+\frac{\zeta_{1}(t)}{\sqrt{2}},$ (40)
$\displaystyle\dot{\rho}_{2}=V\rho_{1}^{2}\sin\theta-\gamma_{2}\rho_{2}\,,$
(41)
$\displaystyle\dot{\theta}={\rho_{1}^{2}-4\rho_{2}^{2}\over\rho_{2}}V\cos\theta+{\sqrt{2}\zeta_{2}(t)\over\rho_{1}}.$
(42)

Here $\zeta_{1}$ and $\zeta_{2}$ are two independent real white noises with
zero mean values $\langle\zeta_{i}(t)\rangle=0$, and the pair correlator
$\langle\zeta_{i}(t_{1})\zeta_{j}(t_{2})\rangle=P\delta_{ij}\delta(t_{1}-t_{2})$.
In the limit of weak noise,
$P|V|^{2}\ll\gamma_{1}^{2}(\gamma_{2}-2\gamma_{1})$, one can apply a linear
approximation near the fixed point. More specifically, we substitute
decomposition
$\rho_{1}(t)=\bar{\rho}_{1}+u(t),\rho_{2}(t)=\bar{\rho}_{2}+v(t),\theta(t)=\bar{\theta}+\phi(t)$
into Eqs. (41)-(42) and keep only the first order terms with respect to $u$,
$v$ and $\phi$. This procedure yields

$\displaystyle\dot{u}=-2V\bar{\rho}_{1}v+{P\over
4\bar{\rho}_{1}}+\frac{\zeta_{1}(t)}{\sqrt{2}},$ (43)
$\displaystyle\dot{v}=2V\bar{\rho}_{1}u-\gamma_{2}v,$ (44)
$\displaystyle\dot{\phi}={(2\gamma_{1}-\gamma_{2})}\phi+{\sqrt{2}\zeta_{2}(t)\over\bar{\rho}_{1}}.$
(45)

From Eqs. (43)-(45) we immediately find the variances

$\displaystyle\langle u^{2}\rangle-\langle
u\rangle^{2}=\frac{\left(2\gamma_{1}+\gamma_{2}\right)P}{8\gamma_{1}\gamma_{2}}\,,$
(46) $\displaystyle\langle v^{2}\rangle-\left\langle
v\right\rangle^{2}=\frac{P}{4\gamma_{2}}\,,$ (47)
$\displaystyle\langle\phi^{2}\rangle-\langle\phi\rangle^{2}=\frac{2V^{2}P}{\gamma_{1}\gamma_{2}\left(\gamma_{2}-2\gamma_{1}\right)}.$
(48)

We see that level of fluctuations in relative phase grows when one approaches
the stability threshold. Note also that far from the threshold, i.e. at
$2\gamma_{1}\ll\gamma_{2}$, the fluctuations of the second harmonic are
suppressed: $\langle v^{2}\rangle/\langle
u^{2}\rangle\approx\gamma_{1}/\gamma_{2}\ll 1$. In this case, the noise of the
first harmonic only weakly influences the conversion into the second one.
However, the conversion is least effective in this limit:
$\bar{\rho}_{2}^{2}/\bar{\rho}_{1}^{2}=\gamma_{1}/2\gamma_{2}\ll 1$.

##  V Conclusion

Our most important finding is the explicit formula (32) for the singular
measure of a direct cascade in the limit of strong noise and weak interaction.
We believe that this is a meaningful advance in non-equilibrium statistics, as
it opens a window to the study of the formation of singular measures in
systems driven far away from equilibrium. We have described also the approach
to this limit and have shown that the total entropy decays and the inter-mode
mutual information grows logarithmically with the Reynolds number. In the
inverse cascade case in this limit, the phase is locked on $-\pi/2$ and the
system exhibits an intermittent dynamics of bursts, which we were able to
describe analytically. Such bursts are perhaps responsible for pulses running
in shell models, which are chains of interacting triplets used in modeling
hydrodynamic incompressible turbulence. It is thus may be interesting to apply
the methods developed here to the popular shell model
$\dot{u}_{i}=u_{i-1}^{2}-u_{i}u_{i+1}$ DN . After some elementary
transformations, this model can be turned into that with the Hamiltonian
${\cal H}=\sum_{i}V_{i}(a_{i}^{2}a_{i+1}^{*}+c.c)$, that is the interacting
chain built of our pairs.

The opposite limit of weak noise and strong interaction is a singular one: the
probability distribution is not close to a Gaussian distribution determined by
$N$ for however small $\chi$, despite occupation numbers being close to
equipartition and the marginal one mode distributions close to quadratic. It
is expressed, in particular, in the nonzero mutual information $I_{12}(\chi)$
at the limit $\chi\to+0$. Figure 1 combines the mutual information data for
both cascades. We see that $I_{12}(\chi)-I_{12}(+0)\propto\chi^{2}$ at
$\chi\ll 1$. We failed to find an analytic solution in this limit either in
the direct or inverse cascade, even though it is likely that the probability
distribution can be expressed in terms of $N$ and ${\cal H}^{2}$, which are
the conserved quantities of the unforced undamped system.

We conclude with suggesting an interesting application of our model to wave
turbulence. In a set of $M+2$ interacting waves, one may consider to model the
interaction of a resonant couple with the other $M$ waves by dissipation and
random forcing. When $M\gg 1$ we can treat forces from all other modes as a
white noise, so that our model (4) applies. In this case, different limits in
$\chi$ correspond to different situations. If we assume an almost continuous
distribution of other modes and estimate from the wave kinetic equation
$\gamma\simeq V^{2}Mn/\omega$ and $P\simeq\gamma n$ ZLF , then
$\chi=\gamma^{3}/PV^{2}\simeq V^{2}Mn/\omega^{2}\ll 1$, which is the original
parameter of nonlinearity assumed to be small. In this case, we come to the
surprising conclusion that a resonant mode within turbulence, when $\Delta
T\simeq T$, has a relative entropy of order unity and independent of $V$. If,
however, we have a set of well-isolated resonant interactions, then it makes
more sense to assume that the interaction with a given mode is $M$ times
smaller so that $\chi$ is large (as $M$ or $\sqrt{M}$), then the relative
entropy is small. Note that in most cases the number of resonant interactions,
is much less than the total number of modes in the system.

The work was supported by the Scientific Excellence Center and Ariane de
Rothschild Women Doctoral Program at WIS, grant 662962 of the Simons
foundation, grant 075-15-2019-1893 by the Russian Ministry of Science, grant
873028 of the EU Horizon 2020 programme, and grants of ISF, BSF and Minerva.
NV was in part supported by NSF grant number DMS-1814619. This work used the
Extreme Science and Engineering Discovery Environment (XSEDE), which is
supported by NSF grant number ACI-1548562, allocation DMS-140028.

## References

  * [1] PD Drummond, KJ McNeil, and DF Walls.  Non-equilibrium transitions in sub/second harmonic generation.  Optica Acta: International Journal of Optics, 27(3):321–335, 1980\.
  * [2] Paul Mandel and Thomas Erneux.  Amplitude self-modulation of intracavity second-harmonic generation.  Optica Acta: International Journal of Optics, 29(1):7–21, 1982\.
  * [3] Matthias W Klein, Christian Enkrich, Martin Wegener, and Stefan Linden.  Second-harmonic generation from magnetic metamaterials.  Science, 313(5786):502–504, 2006.
  * [4] Natalia Vladimirova, Michal Shavit, and Gregory Falkovich.  Fibonacci turbulence.  arXiv preprint arXiv:2101.10418, 2021.
  * [5] Jay Robert Dorfman.  An introduction to chaos in nonequilibrium statistical mechanics.  Number 14. Cambridge University Press, 1999.
  * [6] Michal Shavit and Gregory Falkovich.  Singular measures and information capacity of turbulent cascades.  Physical Review Letters, 125(10):104501, 2020.
  * [7] K Vijay Kumar, Sriram Ramaswamy, and Madan Rao.  Active elastic dimers: Self-propulsion and current reversal on a featureless track.  Physical Review E, 77(2):020102, 2008.
  * [8] Federico Bonetto, Joel L Lebowitz, and Jani Lukkarinen.  Fourier’s law for a harmonic crystal with self-consistent stochastic reservoirs.  Journal of statistical physics, 116(1-4):783–813, 2004.
  * [9] Sergio Ciliberto, Alberto Imparato, Antoine Naert, and Marius Tanase.  Heat flux and entropy produced by thermal fluctuations.  Physical review letters, 110(18):180601, 2013.
  * [10] Gianmaria Falasco, Marco Baiesi, Leo Molinaro, Livia Conti, and Fulvio Baldovin.  Energy repartition for a harmonic chain with local reservoirs.  Physical Review E, 92(2):022129, 2015.
  * [11] Hyun-Myung Chun and Jae Dong Noh.  Hidden entropy production by fast variables.  Physical Review E, 91(5):052128, 2015.
  * [12] Federica Mura, Grzegorz Gradziuk, and Chase P Broedersz.  Nonequilibrium scaling behavior in driven soft biological assemblies.  Physical review letters, 121(3):038002, 2018.
  * [13] Junang Li, Jordan M Horowitz, Todd R Gingrich, and Nikta Fakhri.  Quantifying dissipation using fluctuating currents.  Nature communications, 10(1):1–9, 2019.
  * [14] Paul A Milewski, Esteban G Tabak, and Eric Vanden-Eijnden.  Resonant wave interaction with random forcing and dissipation.  Studies in Applied Mathematics, 108(1):123–144, 2002.
  * [15] RE Lee DeVille, Paul A Milewski, Ricardo J Pignol, Esteban G Tabak, and Eric Vanden-Eijnden.  Nonequilibrium statistics of a reduced model for energy transfer in waves.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 60(3):439–461, 2007.
  * [16] Khanh-Dang Nguyen Thu Lam and Jorge Kurchan.  Stochastic perturbation of integrable systems: a window to weakly chaotic systems.  Journal of Statistical Physics, 156(4):619–646, 2014.
  * [17] VN Desnyansky and EA Novikov.  Modelling of the cascade processes in a turbulent flows.  Prikl. Mat. Mekh, 38(3):507, 1974.
  * [18] Vladimir E Zakharov, Victor S L’vov, and Gregory Falkovich.  Kolmogorov spectra of turbulence I: Wave turbulence.  Springer Science & Business Media, 2012.
  * [19] Natalia Vladimirova, Stanislav Derevyanko, and Gregory Falkovich.  Phase transitions in wave turbulence.  Physical Review E, 85(1):010101, 2012.
  * [20] G Falkovich, I Kolokolov, V Lebedev, and A Migdal.  Instantons and intermittency.  Physical Review E, 54(5):4896, 1996.
  * [21] Gregory Falkovich and Vladimir Lebedev.  Single-point velocity distribution in turbulence.  Physical review letters, 79(21):4159, 1997.
  * [22] AM Obukhov.  Some general properties of equations describing the dynamics of the atmosphere(atmosphere hydrodynamic simulation model for cascade energy transfer in turbulent flow, using euler gyro equations).  Academy of Sciences, USSR, Izvestiya, Atmospheric and Oceanic Physics, 7:471–475, 1971.
  * [23] AJ Bray.  Random walks in logarithmic and power-law potentials, nonuniversal persistence, and vortex dynamics in the two-dimensional xy model.  Physical Review E, 62(1):103, 2000.
  * [24] Artem Ryabov, Ekaterina Berestneva, and Viktor Holubec.  Brownian motion in time-dependent logarithmic potential: Exact results for dynamics and first-passage properties.  The Journal of chemical physics, 143(11):114117, 2015.
  * [25] Martin R Evans and Satya N Majumdar.  Diffusion with stochastic resetting.  Physical review letters, 106(16):160601, 2011.

##  Appendix A Appendix

###  A.1 Hamiltonian evolution

Here we briefly sketch some of the elementary properties of the Hamiltonian
system defined by Eq. (1). From Eq. (1) we find the system of two coupled
complex equations

$\displaystyle\dot{a}_{1}=-i\frac{\partial{\cal H}}{\partial
a_{1}^{*}}=-i\omega a_{1}-2iVa_{1}^{*}a_{2},$ (49)
$\displaystyle\dot{a}_{2}=-i\frac{\partial{\cal H}}{\partial
a_{2}^{*}}=-2i\omega a_{2}-iV^{*}a_{1}^{2}.$ (50)

It is easy to see that in addition to ${\cal H}$, Eqs. (49) and (50) have the
second integral of motion,

$N=\omega|a_{1}|^{2}+2\omega|a_{2}|^{2},$ (51)

and, thus, the system is completely integrable. Indeed, two integrals of
motion allows one to reduce Eqs. (49) and (50) to a single first-order
equation, which we write for $x=2|a_{2}|^{2}/N$ and $t\to tVN^{1/2}$

${dx\over dt}=\pm 2\sqrt{2x(1-x)^{2}-4K^{2}/N^{3}}\ ,$ (52)

where

$K={\cal H}-N=Va_{1}^{*2}a_{2}+V^{*}a_{1}^{2}a_{2}^{*},$ (53)

is also an integral of motion.

Figure 5:  The phase portrait of the integrable Hamiltonian dynamics plotted
for different values of the dimensionless ratio of integrals of motion,
$K^{2}/N^{3}$.

Next, using the Euler representation, $a_{1}=\rho_{1}e^{i\varphi_{1}}$ and
$a_{2}=\rho_{2}e^{i\varphi_{2}}$, one obtains from Eqs. (49) and (50)

$\displaystyle\dot{\rho}_{1}=-2|V|\rho_{1}\rho_{2}\sin\theta\,,$ (54)
$\displaystyle\dot{\rho}_{2}=|V|\rho_{1}^{2}\sin\theta\,,$ (55)
$\displaystyle\dot{\theta}={\rho_{1}^{2}-4\rho_{2}^{2}\over\rho_{2}}\cos\theta,$
(56)

where $\theta=\arg a_{1}^{2}a_{2}^{*}=2\phi_{1}-\phi_{2}$. From Eqs. (52) -
(56) one finds that the dynamical system has two fixed points: 1)
$x=1=2\rho_{2}^{2}/N$, $K=0,$ which means $\rho_{1}=0$ and $\theta=\pm\pi/2$,
and 2) $x=1/3$, $K^{2}/N^{3}=2/27$, which means
$\rho_{1}^{2}=4\rho_{2}^{2}=2N/3$ and $\theta=0$ or $\theta=\pi$. It is easy
to show that the first one is unstable, while the second one is stable. The
small oscillations near the second point are harmonic with the frequency
$4\sqrt{2}/3$. We will see below that this phase portrait explains
qualitatively the statistics of system in the presence of dumping and (small)
random forcing.

Since $K=2V\rho_{1}^{2}\rho_{2}\cos\theta$ is a constant, it cannot change
sign, an, therefore, we have two separated regions in phase space that
correspond to the sign of $\cos\theta$: $\left[-\pi/2,\pi/2\right]$ and
$\left[\pi/2,3\pi/2\right]$. Both separating planes at $\theta=\pi/2,3\pi/2$
are critical points of $\phi$, on which non-linear interaction is zero.

###  A.2 Small-flux limit

This limit can be called alternatively either small-flux limit, because the
mean value of the inter-mode energy flux
$2\langle\rho_{1}^{2}\rho_{2}\sin\theta\rangle$ is much smaller than its
periodic oscillation, or small-noise limit, because most of the time the
evolution is unaffected by pumping and damping. The phase portrait described
in Section A.1 explains qualitatively the statistics of turbulence in this
limit of large wave amplitudes, presented in the upper row of Figure 3.
Indeed, the dynamical system has two fixed points: the unstable one
$\rho_{1}=0$, $\theta=\pm\pi/2$, and the stable one $\rho_{1}^{2}=4\rho$,
$\theta=0,\pi$. In this limit, the system spends most of its time close to one
of the fixed points, randomly switching between them. In the direct cascade,
described in Sect III.1.2, random noise acts on the first mode, so that the
system spends less time around $\rho_{1}$ close to zero, and the probability
has minima at $\theta=\pi/2,-\pi/2$. In the direct cascade, the system spends
more time oscillating around the second fixed point and the probability has
maxima at $\theta=0,\pi$. On the contrary, dissipation acts on $\rho_{1}$ in
the inverse cascade, which keeps system longer around the first fixed point,
and the probability has maxima at $\theta=\pi/2,-\pi/2$.

###  A.3 Large-flux limit for the direct cascade

Here we briefly discuss the derivation of the solution Eq. (36) for the
marginal probability distribution for $\rho_{1}$ in the main text.

Let us plug Eq. (28) into Eq. (17) to get:

$\dot{b}_{1}=-\frac{2|V|^{2}}{\gamma}\frac{|b_{1}|^{2}b_{1}}{1-\frac{4|V|^{2}}{\gamma^{2}}|b_{1}|^{2}}+\frac{\xi(t)}{1-\frac{4|V|^{2}}{\gamma^{2}}|b_{1}|^{2}}.$
(57)

The steady-state probability density of the amplitude $\rho_{1}$ obeys the
Fokker-Planck equation

$\displaystyle\frac{2|V|^{2}}{\gamma}\frac{1}{\rho_{1}}\frac{\partial}{\partial\rho_{1}}\left[\frac{\rho_{1}^{3}}{1-\frac{4|V|^{2}}{\gamma^{2}}\rho_{1}^{2}}{\cal
P}(\rho_{1})\right]+$ (58)
$\displaystyle+\frac{P}{4\rho_{1}}\frac{\partial}{\partial\rho_{1}}\left[\rho_{1}\frac{\partial}{\partial\rho_{1}}\left[\frac{\rho_{1}^{-1}}{(1-\frac{4|V|^{2}}{\gamma^{2}}\rho_{1}^{2})^{2}}{\cal
P}(\rho_{1})\right]\right]=0.$ (59)

Solving this equation one arrives at Eq. (36) in the main text.

###  A.4 Large-flux limit for the inverse cascade

The stationary Fokker-Planck equation on ${\cal P}(\rho_{1},\rho_{2})$ is as
follows

${P\over 4}{\partial^{2}{\cal
P}\over\partial\rho_{2}^{2}}+{\partial\over\partial\rho_{2}}\biggl{[}|V|\rho_{1}^{2}-{P\over
4\rho_{2}}\biggr{]}{\cal
P}+{\partial\over\partial\rho_{1}}\biggl{[}\gamma\rho_{1}-2|V|\rho_{1}\rho_{2}\biggr{]}{\cal
P}=0\ .$ (60)

From Eq. (23) we obtain

$\displaystyle\langle\rho_{2}\rangle=\frac{\gamma}{2|V|},$ (61)
$\displaystyle\langle\rho^{n}_{1}\rangle=\frac{2|V|}{\gamma}\langle\rho_{1}^{n}\rho_{2}\rangle,$
(62)
$\displaystyle\langle\rho_{1}^{2}\rangle=\frac{P}{4|V|}\langle\frac{1}{\rho_{2}}\rangle,$
(63)

Multiplying FPE (60) by $\rho_{2}^{2}$ and integrating over
$d\rho_{1}d\rho_{2}$ yields

$\langle\rho_{1}^{2}\rho_{2}\rangle=\frac{P}{2|V|},$ (64)

and therefore (due to Eqs. (62), (63) and (64))

$\displaystyle\langle\rho_{1}^{2}\rangle=\frac{P}{\gamma},$ (65)
$\displaystyle\langle\frac{1}{\rho_{2}}\rangle=\frac{4|V|}{\gamma}.$ (66)

Also from Eq. (60) we find

$\langle\rho_{2}^{n}\rangle=\frac{4|V|}{(n+2)P}\langle\rho_{1}^{2}\rho_{2}^{n+1}\rangle,\quad\langle\rho_{1}^{2}\rho_{2}^{2}\rangle=\frac{3}{8}\frac{\gamma
P}{|V|^{2}}.$ (67)

It follows from Eqs. (61), (62), (64) and (65) that

$\langle\rho_{1}\rho_{2}\rangle=\langle\rho_{1}\rangle\langle\rho_{2}\rangle,\
\ \
\langle\rho_{1}^{2}\rho_{2}\rangle=\langle\rho_{1}^{2}\rangle\langle\rho_{2}\rangle,$
(68)

which may lead one to hypothesize that in the steady state the random
variables $\rho_{1}$ and $\rho_{2}$ are statistically independent. If such
statistical independence was true, then

${\cal
P}(\rho_{1},\rho_{2})=C\rho_{1}^{-1+\frac{4P|V|^{2}}{\gamma^{3}}}e^{-\frac{2|V|^{2}\rho_{1}^{2}}{\gamma^{2}}}\rho_{2}e^{-\frac{4|V|}{\gamma}\rho_{2}}.$
(69)

However, direct substitution of Eq. (69) into the Fokker-Planck equation (60)
shows that this distribution represents the solution only along two lines:
$\rho_{2}=(\frac{1}{2}+\frac{1}{2\sqrt{2}})\frac{\gamma}{|V|}$ and
$\rho_{1}=\sqrt{\frac{P}{\gamma}}$. Thus, the hypothesis of statistical
independence is not self-consistent.


# Manipulating Predictions over Discrete Inputs in Machine Teaching

Xiaodong Wu  Queen’s University at Kingston, Canada
<EMAIL_ADDRESS>Yufei Han  INRIA, France
<EMAIL_ADDRESS>Hayssam Dahrouj  University of Sharjah, UAE
<EMAIL_ADDRESS>Jianbing Ni  Queen’s University at Kingston, Canada
<EMAIL_ADDRESS>Zhenwen Liang  Notre Dame, America
<EMAIL_ADDRESS>Xiangliang Zhang  Notre Dame, America
<EMAIL_ADDRESS>
###### Abstract

Machine teaching often involves the creation of an optimal (typically minimal)
dataset to help a model (referred to as the ‘student’) achieve specific goals
given by a teacher. While abundant in the continuous domain, the studies on
the effectiveness of machine teaching in the discrete domain are relatively
limited. This paper focuses on machine teaching in the discrete domain,
specifically on manipulating student models’ predictions based on the goals of
teachers via changing the training data efficiently. We formulate this task as
a combinatorial optimization problem and solve it by proposing an iterative
searching algorithm. Our algorithm demonstrates significant numerical merit in
the scenarios where a teacher attempts at correcting erroneous predictions to
improve the student’s models, or maliciously manipulating the model to
misclassify some specific samples to the target class aligned with his
personal profits. Experimental results show that our proposed algorithm can
have superior performance in effectively and efficiently manipulating the
predictions of the model, surpassing conventional baselines.

##  1 Introduction

Machine teaching is a learning paradigm that involves the construction of an
optimal dataset by a teacher model, which enables a student model to learn and
achieve the teacher’s desired objectives. It can be deemed as an inverse
problem to machine learning, which is first proposed by Goldman and Kearns
(1995). Due to its effective application in dataset organization, it has
attracted great interests from the academia. Many previous work (Liu et al.,
2016; Zhu, 2013) study the teaching dimension in the machine teaching problem,
which is the size of the smallest organized dataset for students to learn a
valid model. In addition to the conventional approach of having a student
model mimic the behavior of a teacher model with small datasets, machine
teaching can have an alternative setting where the prediction of a student
model is manipulated by the teacher model.

In the beneficial setting, the teacher leverages the constructed dataset to
guide the student model to behave in a more accurate or robust manner (Han and
Zhang, 2020). For example, some technology companies may seek assistance from
domain experts to help improve their poorly behaved machine learning models.
Teachers in such setting acts like DUTI (Debugging Using Trusted Items). In
another example, Zhang et al. (2018) propose to change training labels as few
as possible by relax the improvement problem into a continuous optimization
problem. Conversely, the malicious setting resembles a targeted poisoning
attack and has recently garnered noticeable attention (Alfeld et al., 2016;
Shafahi et al., 2018). Some prior works (Aghakhani et al., 2021) apply feature
collision attacks to tamper the images so that the predictions of these images
can be manipulated. For both settings, there are already abundant research
works in continuous domain, e.g., images, because the continuous features
(pixels) therein can be modified directly based on the same continuous
gradients, while few attention has been paid to the discrete domain, which is
a critical problem, considering security risks of current powerful nature
language models.

We specifically investigate machine teaching in the context of systems that
rely on discrete inputs in both setting, which are prevalent in machine
learning enabled systems, such as the clinical prediction systems built on
Electronic Health Records (EHR) (Ma et al., 2017). Despite the importance of
this field, it remains relatively unexplored, primarily due to the inherent
challenges associated with discrete inputs. Firstly, machine teaching in the
discrete domain can be formulated as a combinatorial optimization problem,
which has exponentially expanding search space, resulting in the difficulty in
constructing a proper dataset. Secondly, the existing methods in the
continuous domain which use gradient information of the target model to guide
data poisoning cannot be directly applied on discrete training data. This
limitation arises from the fact that the new inputs calculated from the
continuous gradient vectors and original discrete inputs should remain
discrete.

Furthermore, the selection of base samples in continuous domain can be done
randomly, but in discrete domain, a deliberate and meticulous approach is
imperative. Specifically, the samples from the same class in continuous domain
often have similar features, thereby enabling smooth manipulation to construct
an effective teacher dataset. However, applying the same level of perturbation
on different categorical variables can lead to dramatically different
prediction results, as the importance of categorical variables varies.
Therefore, selecting appropriate base samples and determining which features
to perturb become a crucial challenge. A related study by Wang et al. (2020)
propose an optimization-based selection algorithm to solve the feature
selection problem in an evasion attack scenario, where the trained model
remains unchanged. This algorithm tackles the challenge of identifying the
most influential features to perturb for achieving the desired attack
objectives.

Our objective is to construct a compact teacher dataset that, when added to
the training data, can effectively manipulate the predication of the student
model, either correcting erroneous predictions, or leading the model to make
incorrect predictions. The key challenge is to efficiently find a minimized
teacher dataset that achieves the desired manipulation. To address this
challenge, we propose an iterative algorithm called Discrete Machine Teaching
(DMT) that formulates the teacher dataset construction as a combinatorial
optimization problem and tries to solve it with two steps. 1) Teacher dataset
construction. We initially select instances from the training dataset to form
the base of the teacher dataset. For each instance, we determine modifications
to be made on categorical attributes. This involves adding new features,
removing original features, or replacing the original ones with new ones; and
2) student model update.  We combine the teacher dataset with the current
training data and retrain the model. The new model is evaluated to examine if
the manipulation target is reached. If not, we return to step 1) and repeat.
As demonstrated in (Liu et al., 2017), this iterative machine teaching
approach can generate smaller teacher dataset and facilitates faster
convergence. The reason is that this particular approach allows the teacher to
communicate with the student model in multiple rounds. In each round, the
teacher can observe the status of the student and intelligently choose the
most suitable samples to guide the training process towards the desired
target. The key innovation of our DMT lies in the first step. We design
effective score functions for selecting the base samples and features to be
perturbed. Based on the score functions, we achieve successful manipulations
while significantly reducing the searching space. Our code is shared in
Appendix 1.

The main contributions of our work can be summarized as follows.

  * •

We propose to solve the machine teaching problem with discrete categorical
inputs, which is still an under-explored problem.

  * •

We design an efficient method for constructing the teacher dataset, which
estimates the impact of individual samples and their categorical features to
select the most influential samples and features for perturbation. We
incorporate the teacher dataset in an iterative machine teaching framework,
enabling the effective manipulation of student model prediction.

  * •

We evaluate our DMT method on three discrete datasets and observe successful
manipulations on the predictions. Furthermore, our evaluations reveal that our
method outperforms the baseline methods in terms of efficiency.

In the remainder of this paper, we first discuss the related works and their
limitations. Then, we define our problem by providing the problem formulation.
Next, we discuss the details of our methodology. Finally, we evaluate our
method on three discrete datasets and analyze the results.

##  2 Related Work

This section reviews the related work about machine teaching and how to
improve or tamper the predictions in machine teaching.

###  2.1 Machine Teaching

Machine teaching, a new machine learning paradigm, was proposed by Goldman and
Kearns (1995). It aims at creating an optimal dataset from which the student
model learns to meet the goal of the teacher. Subsequently, fundamental
definitions and concepts of machine teaching were further developed in (Zhu,
2015). The current focus of machine teaching study lies in the teaching
dimension, which refers to the size of the minimal dataset required to teach
students to attain the desired objectives (Chen et al., 2018; Kane et al.,
2017). It is desired to have a minimized teaching dimension while achieving
the training target. Previous researches found that introducing interactions
between teacher and student model can effectively reduce the teaching
dimension. For example, Liu et al. (2017) showed how iterative interactions
between the teacher and the student can enhance the process. The iterative
approach significantly reduces the teaching dimension and accelerates
convergence for students. Typically, this method is applied in the continuous
domain, where the optimization problem can be solved with analytical
expressions. However, when used in the discrete domain, the optimization
problem cannot satisfy KKT conditions and generally becomes an NP-hard
problem. Therefore, we propose an efficient iterative machine teaching method
that incorporates an efficient searching strategy to craft an optimal dataset
for discrete data.

###  2.2 Prediction Manipulation

According to the goals of teachers, teachers’ impact on student models can be
classified into two types: prediction improvement and prediction tampering. In
the first case, teacher behaves like DUTI, aiming to improve the
classification prediction results by modifying training items. This debugging
approach has been extensively studied with machine learning techniques
(Cadamuro et al., 2016; Bhadra and Hein, 2015). For example, Han and Zhang
(2020) designed a collaborative machine teaching method that solves the
problem by tuning the training samples with a few trusted items. They optimize
the training samples by utilizing the gradient information from correct
classification, resulting in significant improvement of the student models’
performance.

In the second case, the teacher’s objective is similar to a poisoning attack
where the teacher manipulates the student models to do wrong predictions by
tampering with model’s prediction. Studies in this field are abundant in the
continuous domain. However, in discrete domains, existing works primarily
focus on solving backdoor attack problem, a special type of poisoning attack
(Guo et al., 2022; Salem et al., 2022). The goal of the backdoor attack is to
establish a new relation between a specific target label and some pre-designed
triggers during training. This relationship leads to test samples containing
the pre-designed triggers being classified as the target class. For example,
Schuster et al. (2020) and Yang et al. (2021) explore controlling the word
embeddings within target models. Besides, some other works aim to make their
attacks stealthier and harder to defend against. For instance, Wallace et al.
(2020) propose to using poisoning sentences that do not contain trigger
phrases. This approach achieves high attack success rate and is very difficult
to detect because of the hidden trigger words.

All the above methods are either in the continuous domain or in the discrete
domain but requiring access to tamper both the training and testing data
(backdoor attack). In the general setting, the focus is solely on creating a
perturbed training dataset. Therefore, there is a lack of exploration in
discrete domains where new relation between features and target can be
directly established without relying on triggers. Considering the vast
applications of machine teaching in discrete domain, manipulating student
model’s performance in various NLP tasks is a meaningful but challenging
problem. Our work focuses on generating optimal datasets for efficient machine
teaching in the discrete domain.

##  3 Methodology

In this section, we formulate our problem of teacher dataset construction and
propose the discrete machine teaching methodology to manipulate the
performance of the student models.

###  3.1 Problem Formulation

Let $D_{clean}$ be a clean training dataset without any perturbation, and
contain a set of training data samples $\\{x_{i},y_{i}\\}$, where $y_{i}$ is
the label of $x_{i}$. A classifier $f_{\theta}$ is learned from $D_{clean}$ by
minimizing the objective function
$\mathcal{L}(D_{clean},\theta)=\sum_{(x_{i},y_{i})\in
D_{clean}}l(f_{\theta}(x_{i}),y_{i})$, where $l(.)$ is a typical
classification loss, e.g., cross entropy. The obtained optimal parameter
$\theta_{clean}^{*}$ is expected to make correct predictions in the testing
dataset.

The teacher aims to develop a teacher dataset $D_{teacher}$ such that the
learned $\theta^{*}$ can manipulate prediction on target samples, i.e.,
$f_{\theta^{*}}(x_{j})=\tilde{y}_{j}$, where $\tilde{y}_{j}$ is the target
label that the teacher would like $f_{\theta^{*}}$ to produce with the input
of the target sample $x_{j}$. In addition, the teacher dimension (the size of
$D_{teacher}$) should be minimal. In our discrete setting, $D_{teacher}$ is
constructed by combining the training data $D_{clean}$ with the dataset
$D_{perturbed}$ perturbed from $D_{base}$, which is a selected subset of
$D_{clean}$. Therefore, the overall discrete machine teaching problem is
defined by the objective function as follows,

$\begin{gathered}\mathop{\min}_{D_{perturbed}}\mathcal{L}(D_{target};\mathop{\arg\min}_{\theta}\mathcal{L}(D_{clean}\cup
D_{perturbed};\theta))\\\
s.t.\quad\text{diff}(x_{i},\hat{x}_{i})\leq\epsilon,\;x_{i}\in
D_{base},\hat{x}_{i}\in D_{perturbed},\end{gathered}$ (1)

where $\text{diff}(x_{i},\hat{x}_{i})\leq\epsilon$ establishes the
perturbation budget, $D_{target}$ is the set of targeted samples and labels.
Suppose a discrete instance $x$ is represented by transforming categorical
features into binary variables, i.e., $x\in\\{0,1\\}^{M\times N}$, where $M$
is the number of categorical features and $N$ is the number nominal values
that each feature has. A value 1 in $x$ indicates the presence of the value
within the categorical feature, and 0 means absence. The perturbed $\hat{x}$
can be obtained by changing 1 to 0 in $x$ (_deletion_), or changing 0 to 1 in
$x$ (_insertion_). It is also possible to make several changes and have
_substitution_ in $x$. When $M\times N$ is large, it would be expensive to
take exhaustive search for finding the minimal set of changes that optimizes
Eq. (1).

###  3.2 Teacher’s Capability

The goal of teachers in machine teaching is to shape student model’s behavior,
either by improving their performance in specific tasks, or by stealthily
tampering students’ predictions to benefit teachers. To achieve this goal, the
teachers first should have the necessary computational resources and time
construct their teacher datasets and possession the following knowledge: 1)
The information about the training data used to train the student model. 2)
The architecture of the student model. 3) The details of the target samples
and their expected predictions. Secondly, they should have the access to the
training data containing specific target samples that they aim to influence.
Thirdly, teachers can modify a portion of the training data or insert
generated samples into the dataset, but adhering to a manipulative budget
$\epsilon$. However, they do not have the capacity to modify the test data.
Fourthly, teachers can obtain gradient vectors or middle outputs of the
student model through querying, which guides their modification process. We
next introduce our solution for discrete machine teaching, which includes the
finding of $D_{base}$, and the efficient constructions of $D_{perturbed}$ and
$D_{teacher}$.

###  3.3 DMT Method

Algorithm 1 Discrete Machine Teaching

Input: Original clean dataset $D_{clean}$, target sample $x_{j}$, the
corresponding target class $\tilde{y}_{j}$, a student model $f_{\theta}(.)$,
the max number of iteration $T$

Output: Manipulated model $f_{\theta^{*}}(.)$

1:
$\theta^{*}=\operatorname*{argmin}\limits_{\theta}\mathcal{L}(D_{clean};\theta)$;

2: $D_{0}=D_{clean}$;

3: for $t=1$ to $T$ do

4: $D_{base}=k$-nn$(x_{j},D_{t-1})$;

5: $D_{perturbed}=GGGM(D_{base},x_{j},\theta^{*})$;

6: $D_{t}=Combine(D_{t-1},D_{perturbed})$;

7:
$\theta^{*}=\operatorname*{argmin}\limits_{\theta}\mathcal{L}(D_{t};\theta)$;

8: if $f_{\theta^{*}}(x_{j})=\tilde{y}_{j}$ then

9: break

10: end if

11: end for

The proposed DMT method is presented in Algorithm 1. It is an iterative
algorithm including basic steps of 1) teacher dataset construction (line 4-6),
and 2) student model update (line 7). The new model is tested if the
manipulation target is reached (line 8). If not, go to step 1) and repeat
until hitting the maximum number of iterations $T$. The construction of the
teacher dataset consists of three steps below.

Selection of $D_{base}$. To minimize $D_{teacher}$, $D_{base}$ should consist
of the most influential base samples that have the capability to alter the
prediction of target sample $x_{j}$. Intuitively, the samples in $D_{base}$
should surround $x_{j}$. We select base samples by determining the $k$-nearest
neighbors of the target sample $x_{j}$ in the training dataset. Since the data
have discrete categorical features, we use the Jaccard distance for the
nearest neighbor identifying. Initially, $D_{base}$ is selected from the clean
dataset $D_{0}=D_{clean}$. When the teacher dataset $D_{t}$ is updated, the
selection of $D_{base}$ in next round would have access to the updated
$D_{t}$, including the perturbed dataset $D_{perturbed}$ at step $t$. The
selected $k$ samples would be used to create $D_{perturbed}$ based on Gradient
Guided Greedy Method (GGGM), which is introduced in the next section. In the
iterative teaching process, $k$ is also named as _step size_. Our evaluation
results demonstrate the impact of $k$ on the teaching performance.

Construction of $D_{perturbed}$. The construction of $D_{perturbed}$ is to
select the most powerful features in each sample of $D_{base}$. The details
are shown in Algorithm 2. For each sample in $D_{base}$, we use gradients to
first filter out the candidate changes that have the largest absolute gradient
values (line 5-6), as the gradient magnitude of $f_{\theta^{*}}(x)$ over $x$
implies the contribution to prediction (Lei et al., 2018). Unlike what was
done in continuous domain, these gradients cannot be used directly to
construct $D_{perturbed}$, but can serve as guiding which feature to be
selected and modified. Then, we evaluate the impact of these candidate changes
by designing a function $g$ and choosing those with strongest potential to
manipulate the prediction of $x_{j}$ (line 7). The perturbations made to
obtain $\hat{x}$ is controlled within the budget $\epsilon$.

To evaluate the impact of modifying a subset of candidate features, we design
two score functions $g(\cdot)$. The first one is based on the Euclidean
distance function $g_{dist}$, which measures the difference of predictions
made on the target sample $x_{j}$ and the perturbed sample $\hat{x}$,

$\begin{gathered}g_{dist}(\hat{x},x_{j},\theta)=\|f_{\theta^{*}}(x_{j})-f_{\theta^{*}}(\hat{x})\|_{2}.\end{gathered}$
(2)

A perturbation is preferred if the perturbed $\hat{x}$ has a closer prediction
score compared with the score of the target sample $x_{j}$.

The second method is the alignment function $g_{align}$ inspired from Geiping
et al. (2020):

$\small\begin{gathered}g_{align}(\hat{x},x,x_{j},\theta^{*})=-\lambda*align(\hat{x},x_{j},\theta^{*})\\\
-(1-\lambda)*align(\hat{x},x,\theta^{*}),\end{gathered}$ (3)

where $\lambda$ is a trade-off parameter that combines two alignment scores:
one measures how a perturbed sample $\hat{x}$ is aligned with the targeted
sample $x_{j}$; the other one is measuring the alignment between the perturbed
sample $\hat{x}$ and the clean sample $x$. The alignment score is defined as
the cosine distance of the gradient vectors of two given samples, which is
defined as:

$\small
align(x_{a},x_{b},\theta)=\frac{\left<\bigtriangledown_{\theta}l(f_{\theta}(x_{a}),y_{a}),\bigtriangledown_{\theta}l(f_{\theta}(x_{b}),\tilde{y}_{b})\right>}{\|\bigtriangledown_{\theta}l(f_{\theta}(x_{a}),y_{a})\|\|\bigtriangledown_{\theta}l(f_{\theta}(x_{b}),\tilde{y}_{b})\|},$
(4)

The idea behind this design is to strike a balance between achieving the
manipulation goal (via comparing $\hat{x}$ and $x_{j}$) and controlling the
modification extend (via comparing $\hat{x}$ and clean $x$).

Algorithm 2 Gradient Guided Greedy Method (GGGM)

Input: base dataset $D_{base}$, a trained student model $f_{\theta^{*}}(.)$,
the target sample $x_{j}$, the manipulative budget $\epsilon$, the candidate
set size $q$

Output: perturbed dataset $D_{perturbed}$

1: $D_{perturbed}=\emptyset$

2: for $x$ in $D_{base}$ do

3: $\hat{x}^{1}=x$

4: for $p=1$ to $\epsilon/q$ do

5: $r=\bigtriangledown_{x}f_{\theta^{*}}(\hat{x}^{p})$

6: $s=\\{(m,n)\\}_{|r_{mn}|\in\;\text{top-}q\text{-values}(|r_{.,.}|)}$

7: $S=\operatorname*{argmin}\limits_{l\subset
s}g($modify$(\hat{x}^{p},l),x,x_{j},\theta^{*})$

8: $\hat{x}^{p+1}$ = modify $(\hat{x}^{p},S)$

9: end for

10: $D_{perturbed}\leftarrow D_{perturbed}\cup\hat{x}^{p}$

11: end for

Combination of $D_{perturbed}$ and $D_{base}$. In Algorithm 1 line 6, after
obtaining the new dataset $D_{perturbed}$, a strategy should be applied to
combine this dataset with the base dataset $D_{t-1}$ from the last iteration.
Here, we propose to incrementally add new generated fake data $D_{perturbed}$
to the base dataset $D_{t-1}$. Consequently, the size of the modified dataset
$D_{t}$ increases by a step size $k$, i.e., the cardinality of
$D_{perturbed}$, in each iteration. Due to the nature of this strategy,
incremental addition allows the effect of the manipulation dataset in each
iteration gradually manipulates the prediction of the target sample toward the
target class step by step. Furthermore, by keeping the step size small, the
distribution of $D_{t-1}$ would not be changed dramatically, enabling the
model to predict other test samples with a reasonable accuracy. An alternative
strategy is to _replace_ the base samples in $D_{t-1}$ that correspond to the
samples in $D_{perturbed}$. However, as this strategy creates larger changes
in $D_{t-1}$, the distribution of $D_{t-1}$ would be significantly affected.
Further comparisons between these two strategies, i.e., _addition_ and
_replacement_ , are provided in Appendix 4.

In addition, when there are multiple samples in $D_{target}$, and the teacher
aims to manipulate a group of targeted samples, we calculate the average
distance regarding all target samples in the step of base sample selection
(line 4 of Algorithm 1). Similarly, for line 7 of Algorithm 2, $g_{dist}$ is
calculated by taking the average over all target samples. When using the
alignment function in Eq.(3), we first average the gradients across all target
samples, and then use this averaged gradient to calculate the alignment score.

Figure 1: The CSR of Different Methods on the Performance Improvement Task
When Varying the Percentages of Samples Allowed to Change.

Figure 2:  The CSR of Different Methods in Prediction Tampering Tasks with
Varying Sample Percentage in MALWARE, EHR, and IPS.

##  4 Experiments

In this section, we present the dataset used to evaluate the performance of
the DMT method and demonstrate its effectiveness on both prediction
improvement and tempering tasks.

###  4.1 Dataset

We use three discrete datasets from different application domains to evaluate
the performance of our proposed DMT method, which are EHR, MALWARE and IPS.
The details of them are introduced in Appendix 2.

###  4.2 Implementation Details

The implementation details of our experimental results are shown in Appendix
3.

###  4.3 Evaluation Tasks and Metrics

We evaluate the performance of DMT from two perspectives. One is to use DMT to
improve the performance of a student model, i.e., correcting the wrongly
classified samples. The second perspective is to allow the teacher to tamper
student model’s performance maliciously and share the same problem setting of
poisoning attacks. It is worth noting that the performance improvement task
requires the successful manipulation of the entire group of samples, rather
than just one sample, as in the performance tampering tasks. In other words,
the break condition in Algorithm 1 line 8 requires that all target samples in
one group should be correctly classified. Therefore, the performance
improvement tasks are more challenging than performance tampering tasks.

In our comparison of different methods, we focus on both computational
efficiency and performance. Regarding the efficiency, we evaluate the average
time consumption for each iteration, the average number of iterations required
and the average percentage of samples needed over the entire dataset to
achieve successful manipulation. Furthermore, we examine the change success
rate (CSR), i.e., the number of correct classification of target samples as
their desired class across all target samples. CSRs of different methods are
compared under the constraint of the same maximum percentage of samples
allowed to change. Finally, we consider the average decision logit of the last
5 training epochs as the output target logit to decide whether the target
sample can be classified as the target class.

###  4.4 Baseline Methods

As there is no prior work of discrete machine teaching, we take a variant of
our DMT as one baseline. There are two other models designed for the
continuous domain. We adapt them to manipulate the discrete categorical
features for performance improvement and prediction tampering tasks.

  * •

At once is a method based on our proposed method but with only one iteration.
It can be considered as our method with a big step size $k$. It only inserts
the generated teacher dataset into the clean dataset so as to obtain the
results at once.

  * •

Frogs Shafahi et al. (2018) is a feature collision attack method, which
proposed to move the base samples toward the target sample in the embedding
space, rather than make the perturbation in the original feature space.

  * •

Witch’s brew Geiping et al. (2020) formulated the poisoning attack as a
gradient matching problem and tackled the problem by aligning the target and
poisoning gradients in the same direction.

Apart from the first baseline, the other two methods were proposed to solve
the poisoning attack problem in the continuous domain. Since the input data in
the discrete domain are preprocessed as one hot vectors, we apply a rounding
process to the samples output from these algorithms so that they can be
converted into the discrete domain. For all baseline methods, we test their
performances using different fixed numbers of manipulated samples to make a
comparison.

###  4.5 Prediction Improvement Results

TABLE I: Efficiency Evaluation with Different Score Function on MALWARE and IPS Dataset for Improving Predictions Dataset | Method | Step Size | Iteration | Time(min) | Sample Percent
---|---|---|---|---|---
MALWARE | DMT-Dist | k=10 | 4.7 | 0.7 | 0.24%
k=20 | 4.0 | 0.8 | 0.36%
DMT-Align | k=10 | 2.6 | 1.1 |  0.13%
k=20 | 2.3 | 1.5 | 0.23%
IPS | DMT-Dist | k=40 | 7.2 | 0.4 | 7.58%
k=50 | 7.2 | 0.5 | 9.47%
DMT-Align | k=40 | 6.2 | 3.5 |  6.53%
k=50 | 5.4 | 4.6 | 7.11%

The goal of a teacher is to correct the predictions of a group of wrongly
classified samples. We report the results of our methods and baselines on
MALWARE and IPS. The efficiency evaluation results are given in Table.I and
results of performance improvement are shown in Figure 1.

In Figure 1, the differences among evaluated methods may not be prominent.
However, we can still observe some trends. Notably, the smaller step size
leads to better improvement performances. Our proposed method using alignment
score function outperforms the one using distance function. Besides, compared
to baseline methods, our best setting, i.e., alignment function with small
step size can achieve the highest CSR under the same sample percentage. As the
sample percentage increases, the baseline methods exhibit similar improvement
performance to the DMT-dist method, but the DMT-dist can eventually achieve
100% CSR with a reasonable sample percentage limit.

In terms of efficiency, from Table.I, we can find that a larger step size
results in fewer average iteration and required samples, but it also consumes
more computing time. Compared with the distance method, the alignment score
based method can reduce the number of iterations and generated samples, but it
requires more time to achieve this reduction.

TABLE II: Efficiency Evaluation with Different Score Function on MALWARE, IPS and EHR in Prediction Tampering Tasks DataSet | Method | Step Size | Iteration | Time(min) | Sample Percent
---|---|---|---|---|---
MALWARE | DMT-Dist | k=100 | 14.2 | 8.8 | 7.12%
k=200 | 11.4 | 10.2 | 11.3%
k=400 | 11.1 | 11.0 | 22.28%
DMT-Align | k=100 | 11.4 | 13.7 |  5.72%
k=200 | 9.8 | 23.9 | 9.77%
k=400 | 8.1 | 44.6 | 16.27%
IPS | DMT-Dist | k=40 | 9.8 | 1.2 | 8.47%
k=50 | 8.0 | 1.6 | 9.27%
k=70 | 6.2 | 2.1 | 10.50%
DMT-Align | k=40 | 9.1 | 3.1 |  8.15%
k=50 | 7.3 | 4.1 | 8.43%
k=70 | 5.7 | 5.7 | 9.68%
EHR | DMT-Dist | k=10 | 9.6 | 24.0 | 1.51%
k=20 | 7.1 | 26.9 | 2.21%
k=40 | 4.6 | 32.7 | 2.88%
DMT-Align | k=10 | 7.8 | 26.6 |  1.21%
k=20 | 6.5 | 30.5 | 2.04%
k=40 | 4.1 | 39.8 | 2.58%

###  4.6 Prediction Tampering Results

We show the prediction tampering results of our machine teaching method in
Table II. The teacher plays a role of malicious attackers on tampering the
prediction results on MALWARE, EHR and IPS dataset. The performance of
attacking is shown in Figure 2. We evaluate our DMT using different score
functions: distance and alignment score. For each DMT with a specific score
function, the step size $k$ is set to 3 different values. The curves in Figure
2 illustrate the increase of CSR when more samples are introduced into
$D_{teacher}$ in the iterative teaching process. All curves are compared at
the same sample percentage. A small $k$ leads to a slow increase in
$D_{teacher}$, while a large $k$ results in a faster increase. The baseline
method _at once_ takes the entire set $D_{perturbed}$ at a given sample
percentage into $D_{teacher}$ by one step, without running iteratively.

The comparison results reveal several important findings. First, when the step
size $k$ is smaller, a high CSR can be reached with a smaller percentage of
samples. This indicates that a smaller step size can lead to more powerful
manipulation. This is because a small $k$ value can accumulate the most
powerful samples after enough iterations through interactive k-NN selection.
It is also observed that our method outperforms all baseline methods.
Moreover, comparing our DMT using distances and alignment score functions with
the same step size, the alignment score function can result in greater CSRs.
This aligns with our prediction that poisoned samples selected based on
alignment score, using the gradients from the model are more powerful than
ones selected based on the distance score. The alignment score function is
more effective since it exploits the relation between manipulation samples and
target samples.

Second, with the increase of step size, there is a noticeable rise in time
consumption and the percentage of fake samples required to achieve successful
manipulation. It indicates that no matter what the feature selection strategy
is, a larger step size would lead to higher computational cost. Comparing the
distance based method and the alignment score method, the latter can save more
on sample perturbation and teaching iterations but consuming more time when
crafting these samples. This trade-off between efficiency and performance
should be considered when choosing a specific method for discrete machine
teaching.

##  5 Conclusion

In this paper, we propose an iterative discrete machine teaching algorithm,
involving teacher dataset construction and student model update. This is the
first work investigating machine teaching problem in discrete domain without
access to test data. We exploit score functions to select the most influential
base samples from the dataset and then choose the most important features to
change. Experimental results on three datasets demonstrate the effectiveness
of our DMT in both performance improvement and performance tampering tasks,
achieving strong and efficient manipulation with up to $100\%$ change success
rate and minimal time and samples consumption. The results in this paper
represent one step forward in the area of machine teaching on discrete data,
and open up exciting avenues for future research, e.g., to investigate change
success guarantees for given target samples.

## References

  * Agarap (2018) Abien Fred Agarap.  Deep learning using rectified linear units (relu).  _arXiv preprint arXiv:1803.08375_ , 2018.
  * Aghakhani et al. (2021) Hojjat Aghakhani, Dongyu Meng, Yu-Xiang Wang, Christopher Kruegel, and Giovanni Vigna.  Bullseye polytope: A scalable clean-label poisoning attack with improved transferability.  In _2021 IEEE European Symposium on Security and Privacy (EuroS &P)_, pages 159–178. IEEE, 2021.
  * Akgun et al. (2012) Baris Akgun, Maya Cakmak, Jae Wook Yoo, and Andrea Lockerd Thomaz.  Trajectories and keyframes for kinesthetic teaching: A human-robot interaction perspective.  In _Proceedings of the seventh annual ACM/IEEE international conference on Human-Robot Interaction_ , pages 391–398, 2012.
  * Alfeld et al. (2016) Scott Alfeld, Xiaojin Zhu, and Paul Barford.  Data poisoning attacks against autoregressive models.  In _AAAI_ , volume 30, 2016.
  * Bao et al. (2021) Hongyan Bao, Yufei Han, Yujun Zhou, Yun Shen, and Xiangliang Zhang.  Towards understanding the robustness against evasion attack on categorical data.  In _ICLR_ , 2021.
  * Bhadra and Hein (2015) Sahely Bhadra and Matthias Hein.  Correction of noisy labels via mutual consistency check.  _Neurocomputing_ , 160:34–52, 2015.
  * Cadamuro et al. (2016) Gabriel Cadamuro, Ran Gilad-Bachrach, and Xiaojin Zhu.  Debugging machine learning models.  In _ICML Workshop on Reliable Machine Learning in the Wild_ , volume 103, 2016.
  * Chan et al. (2020) Alvin Chan, Yi Tay, Yew-Soon Ong, and Aston Zhang.  Poison attacks against text datasets with conditional adversarially regularized autoencoder.  _arXiv preprint arXiv:2010.02684_ , 2020.
  * Chen et al. (2018) Yuxin Chen, Adish Singla, Oisin Mac Aodha, Pietro Perona, and Yisong Yue.  Understanding the role of adaptivity in machine teaching: The case of version space learners.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Cheng et al. (2016) Yu Cheng, Fei Wang, Ping Zhang, and Jianying Hu.  Risk prediction with electronic health records: A deep learning approach.  In _Proceedings of the 2016 SIAM international conference on data mining_ , pages 432–440. SIAM, 2016.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * Doliwa et al. (2014) Thorsten Doliwa, Gaojian Fan, Hans Ulrich Simon, and Sandra Zilles.  Recursive teaching dimension, vc-dimension and sample compression.  _The Journal of Machine Learning Research_ , 15(1):3107–3131, 2014.
  * Fang et al. (2020) Minghong Fang, Neil Zhenqiang Gong, and Jia Liu.  Influence function based data poisoning attacks to top-n recommender systems.  In _Proceedings of The Web Conference 2020_ , pages 3019–3025, 2020.
  * Geiping et al. (2020) Jonas Geiping, Liam Fowl, W Ronny Huang, Wojciech Czaja, Gavin Taylor, Michael Moeller, and Tom Goldstein.  Witches’ brew: Industrial scale data poisoning via gradient matching.  _arXiv preprint arXiv:2009.02276_ , 2020.
  * Goldman and Kearns (1995) Sally A Goldman and Michael J Kearns.  On the complexity of teaching.  _Journal of Computer and System Sciences_ , 50(1):20–31, 1995.
  * Goodfellow et al. (2014) Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  _arXiv preprint arXiv:1412.6572_ , 2014.
  * Guo and Liu (2020) Junfeng Guo and Cong Liu.  Practical poisoning attacks on neural networks.  In _ECCV_ , pages 142–158. Springer, 2020.
  * Guo et al. (2022) Wei Guo, Benedetta Tondi, and Mauro Barni.  An overview of backdoor attacks against deep neural networks and possible defences.  _IEEE Open Journal of Signal Processing_ , 2022.
  * Han and Zhang (2020) Yufei Han and Xiangliang Zhang.  Robust federated learning via collaborative machine teaching.  In _AAAI_ , volume 34, pages 4075–4082, 2020.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Jagielski et al. (2021) Matthew Jagielski, Giorgio Severi, Niklas Pousette Harger, and Alina Oprea.  Subpopulation data poisoning attacks.  In _Proceedings of the 2021 ACM SIGSAC Conference on Computer and Communications Security_ , pages 3104–3122, 2021.
  * Kane et al. (2017) Daniel M Kane, Shachar Lovett, Shay Moran, and Jiapeng Zhang.  Active classification with comparison queries.  In _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 355–366. IEEE, 2017.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Koh and Liang (2017) Pang Wei Koh and Percy Liang.  Understanding black-box predictions via influence functions.  In _ICML_ , pages 1885–1894. PMLR, 2017.
  * Koh et al. (2022) Pang Wei Koh, Jacob Steinhardt, and Percy Liang.  Stronger data poisoning attacks break data sanitization defenses.  _Machine Learning_ , 111(1):1–47, 2022.
  * Kouliaridis and Kambourakis (2021) Vasileios Kouliaridis and Georgios Kambourakis.  A comprehensive survey on machine learning techniques for android malware detection.  _Information_ , 12(5):185, 2021.
  * Kurita et al. (2020) Keita Kurita, Paul Michel, and Graham Neubig.  Weight poisoning attacks on pre-trained models.  _arXiv preprint arXiv:2004.06660_ , 2020.
  * Lei et al. (2018) Qi Lei, Lingfei Wu, Pin-Yu Chen, Alexandros G Dimakis, Inderjit S Dhillon, and Michael Witbrock.  Discrete attacks and submodular optimization with applications to text classification.  _arXiv preprint arXiv:1812.00151_ , 2018.
  * Liu et al. (2016) Ji Liu, Xiaojin Zhu, and Hrag Ohannessian.  The teaching dimension of linear learners.  In _International Conference on Machine Learning_ , pages 117–126. PMLR, 2016.
  * Liu et al. (2017) Weiyang Liu, Bo Dai, Ahmad Humayun, Charlene Tay, Chen Yu, Linda B Smith, James M Rehg, and Le Song.  Iterative machine teaching.  In _ICML_ , pages 2149–2158. PMLR, 2017.
  * Ma et al. (2017) Fenglong Ma, Radha Chitta, Jing Zhou, Quanzeng You, Tong Sun, and Jing Gao.  Dipole: Diagnosis prediction in healthcare via attention-based bidirectional recurrent neural networks.  In _KDD_ , KDD ’17, page 1903–1911, 2017.
  * Mei and Zhu (2015) Shike Mei and Xiaojin Zhu.  Using machine teaching to identify optimal training-set attacks on machine learners.  In _AAAI_ , 2015.
  * Muñoz-González et al. (2019) Luis Muñoz-González, Bjarne Pfitzner, Matteo Russo, Javier Carnerero-Cano, and Emil C Lupu.  Poisoning attacks with generative adversarial nets.  _arXiv preprint arXiv:1906.07773_ , 2019.
  * Salem et al. (2022) Ahmed Salem, Rui Wen, Michael Backes, Shiqing Ma, and Yang Zhang.  Dynamic backdoor attacks against machine learning models.  In _2022 IEEE 7th European Symposium on Security and Privacy (EuroS &P)_, pages 703–718. IEEE, 2022.
  * Schuster et al. (2020) Roei Schuster, Tal Schuster, Yoav Meri, and Vitaly Shmatikov.  Humpty dumpty: Controlling word meanings via corpus poisoning.  In _2020 IEEE Symposium on Security and Privacy (SP)_ , pages 1295–1313. IEEE, 2020.
  * Shafahi et al. (2018) Ali Shafahi, W Ronny Huang, Mahyar Najibi, Octavian Suciu, Christoph Studer, Tudor Dumitras, and Tom Goldstein.  Poison frogs! targeted clean-label poisoning attacks on neural networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Simard et al. (2017) Patrice Y Simard, Saleema Amershi, David M Chickering, Alicia Edelman Pelton, Soroush Ghorashi, Christopher Meek, Gonzalo Ramos, Jina Suh, Johan Verwey, Mo Wang, et al.  Machine teaching: A new paradigm for building machine learning systems.  _arXiv preprint arXiv:1707.06742_ , 2017.
  * Wallace et al. (2020) Eric Wallace, Tony Z Zhao, Shi Feng, and Sameer Singh.  Concealed data poisoning attacks on nlp models.  _arXiv preprint arXiv:2010.12563_ , 2020.
  * Wang et al. (2021) Jun Wang, Chang Xu, Francisco Guzmán, Ahmed El-Kishky, Yuqing Tang, Benjamin IP Rubinstein, and Trevor Cohn.  Putting words into the system’s mouth: A targeted attack on neural machine translation using monolingual data poisoning.  _arXiv preprint arXiv:2107.05243_ , 2021.
  * Wang et al. (2020) Yutong Wang, Yufei Han, Hongyan Bao, Yun Shen, Fenglong Ma, Jin Li, and Xiangliang Zhang.  Attackability characterization of adversarial evasion attack on discrete data.  In _KDD_ , pages 1415–1425, 2020.
  * Yang et al. (2017) Chaofei Yang, Qing Wu, Hai Li, and Yiran Chen.  Generative poisoning attack method against neural networks.  _arXiv preprint arXiv:1703.01340_ , 2017.
  * Yang et al. (2021) Wenkai Yang, Lei Li, Zhiyuan Zhang, Xuancheng Ren, Xu Sun, and Bin He.  Be careful about poisoned word embeddings: Exploring the vulnerability of the embedding layers in nlp models.  _arXiv preprint arXiv:2103.15543_ , 2021.
  * Yang et al. (2019) Zhilin Yang, Zihang Dai, Yiming Yang, Jaime Carbonell, Russ R Salakhutdinov, and Quoc V Le.  Xlnet: Generalized autoregressive pretraining for language understanding.  _Advances in neural information processing systems_ , 32, 2019.
  * Zhang et al. (2018) Xuezhou Zhang, Xiaojin Zhu, and Stephen Wright.  Training set debugging using trusted items.  In _AAAI_ , volume 32, 2018.
  * Zhao et al. (2017) Zhengli Zhao, Dheeru Dua, and Sameer Singh.  Generating natural adversarial examples.  _arXiv preprint arXiv:1710.11342_ , 2017.
  * Zhu (2013) Jerry Zhu.  Machine teaching for bayesian learners in the exponential family.  _Advances in Neural Information Processing Systems_ , 26, 2013.
  * Zhu (2015) Xiaojin Zhu.  Machine teaching: An inverse problem to machine learning and an approach toward optimal education.  In _AAAI_ , volume 29, 2015.
  * Zilles et al. (2011) Sandra Zilles, Steffen Lange, Robert Holte, Martin Zinkevich, and Nicolò Cesa-Bianchi.  Models of cooperative teaching and learning.  _Journal of Machine Learning Research_ , 12(2), 2011.