# Anticipating critical transitions in multi-dimensional systems driven by
time- and state-dependent noise

Andreas Morr<EMAIL_ADDRESS>Earth System Modelling, School of
Engineering and Design, Technical University of Munich, Germany  Research
Domain IV – Complexity Science, Potsdam Institute for Climate Impact Research,
Germany  Keno Riechers  Earth System Modelling, School of Engineering and
Design, Technical University of Munich, Germany  Research Domain IV –
Complexity Science, Potsdam Institute for Climate Impact Research, Germany
Leonardo Rydin Gorjão  Faculty of Science and Technology, Norwegian University
of Life Sciences, Norway  Niklas Boers  Earth System Modelling, School of
Engineering and Design, Technical University of Munich, Germany  Research
Domain IV – Complexity Science, Potsdam Institute for Climate Impact Research,
Germany  Global Systems Institute and Department for Mathematics, University
of Exeter, UK

###### Abstract

The anticipation of bifurcation-induced transitions in dynamical systems has
gained relevance in various fields of the natural, social, and economic
sciences. When approaching a co-dimension 1 bifurcation, the feedbacks that
stabilise the initial state weaken and eventually vanish; a process referred
to as critical slowing down (CSD). This motivates the use of variance and
lag-1 autocorrelation as indicators of CSD. Both indicators rely on
linearising the system’s restoring rate. Additionally, the use of variance is
limited to time- and state-independent driving noise, strongly constraining
the generality of CSD. Here, we propose a data-driven approach based on
deriving a Langevin equation to detect local stability changes and anticipate
bifurcation-induced transitions in systems with generally time- and state-
dependent noise. Our approach substantially generalizes the conditions
underlying existing early warning indicators, which we showcase in different
examples. Changes in deterministic dynamics can be clearly discriminated from
changes in the driving noise. This reduces the risk of false and missed alarms
of conventional CSD indicators significantly in settings with time-dependent
or multiplicative noise. In multi-dimensional systems, our method can greatly
advance the understanding of the coupling between system components and can
avoid risks of missing CSD due to dimension reduction, which existing
approaches suffer from.

††preprint: APS/123-QED

##  I Introduction

A mechanistic understanding of complex high-dimensional physical systems is
essential for assessing the risk of abrupt regime shifts, for example in
ecological, climatic, social, or financial systems. Such shifts may occur when
critical forcing thresholds, which correspond to underlying bifurcation
points, are crossed [1, 2, 3, 4]. Reducing complex systems to a low-
dimensional summary observable $\mathbf{X}_{t}$ has leveraged impressive
modelling capabilities [5, 6, 7, 8]. This is particularly important because
observations are typically available in the form of multivariate time series
of just a few dimensions.

Commonly, the dynamics of the summary observable
$\mathbf{X}_{t}\in\mathbb{R}^{n}$ is approximately separated into a
deterministic component $A(\mathbf{X}_{t},t)\mathrm{d}t$ and a stochastic
component $B(\mathbf{X}_{t},t)\mathrm{d}\mathbf{W}_{t}$ that represents the
action of the omitted dimensions. This results in the Langevin equation

$\mathrm{d}\mathbf{X}_{t}=A(\mathbf{X}_{t},t)\mathrm{d}t+B(\mathbf{X}_{t},t)\mathrm{d}\mathbf{W}_{t}.$
(1)

Even though, in principle, the stochastic component can take more complicated
forms, we restrict ourselves to the case where $\mathbf{W}$ is an uncorrelated
Wiener process supported on the filtered probability space
$(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\in\mathbb{R}_{+}},\mathbb{P})$ and
refer to existing extensions to the case of correlated noise [9, 10].

This framework facilitates a mathematical description of abrupt regime shifts
in terms of dynamic bifurcations in low-dimensional dynamical systems [2, 11].
Prior to the transition, the deterministic drift $A(\mathbf{X}_{t},t)$
embodies negative feedback mechanisms keeping the system in a stable
equilibrium [12, 13, 14]. The explicit time dependence of
$A(\mathbf{X}_{t},t)$ reflects the changing forcing levels that act on the
system from the outside and alter the deterministic, coarse-grained dynamics.
At the bifurcation point, i.e. at the critical level of forcing, the currently
occupied equilibrium state is annihilated and the system abruptly transitions
to another stable state.

For co-dimension 1 bifurcations, it is well known that a weakening of the
negative feedback precedes an eventual abrupt transition [2, 15]. This will,
heuristically speaking, result in a weaker and slower response to the pseudo-
random perturbations stemming from the unresolved dynamics. This phenomenon is
referred to as critical slowing down (CSD), and it manifests in an increase of
the statistical quantities of variance and lag-1 autocorrelation (AC(1)) of
the observable in the components exhibiting stability loss [16, 17, 2, 18].
These two quantities are therefore often employed to anticipate bifurcation-
induced abrupt transitions, and their simultaneous increase has been suggested
as an early warning signal (EWS) [2, 19, 3]. Mathematically, CSD can be
described by approximating the negative feedback around a stable equilibrium
as a linear restoring rate. Denoting the time-dependent equilibrium state of a
one-dimensional observable $X_{t}$ by $x^{*}(t)$, we arrive at the
Ornstein–Uhlenbeck model [20]

$\mathrm{d}X_{t}=-\lambda(t)(X_{t}-x^{*}(t))\mathrm{d}t+\sigma\mathrm{d}W_{t}.$
(2)

The linearised negative feedback $\lambda(t)$ weakens during CSD while the
noise-coupling strength $\sigma$ is assumed to remain constant. This results
in the following expressions for variance and AC(1),

$\displaystyle\mathrm{Var}\left[X\right]$
$\displaystyle=\frac{\sigma^{2}}{2\lambda}\xrightarrow{\lambda\rightarrow
0}\infty$ (3) $\displaystyle\mathrm{AC}_{X}(1)$
$\displaystyle=\exp(-\lambda\Delta t)\xrightarrow{\lambda\rightarrow 0}1,$

where $\Delta t>0$ is the sampling time step of the data. Detection of CSD is
usually preceded by a reduction of the system to a one-dimensional observable,
either by leveraging physical understanding or employing principle component
analysis [16, 21] in order to identify a linear combination of components
which may be experiencing stability loss. In such a reduction, crucial
information about system stability may be lost. The method presented herein is
applicable directly to data from higher dimensional systems (or to
multivariate data) and thus avoids this preprocessing step.

We will nevertheless, for illustrative purposes, first treat the problem of
estimating local system stability in the one-dimensional dynamical system
denoted by

$\mathrm{d}X_{t}=a(X_{t},t)\mathrm{d}t+b(X_{t},t)\mathrm{d}W_{t},$ (4)

such that the linearised negative feedback takes the form

$\lambda(t):=-\partial_{x}a(x^{*}(t),t).$ (5)

The extension of the discussed estimation methods to the general, higher-
dimensional setting of Eq. (1) is discussed in the Methods section.

Time- or state-dependent driving noise can lead to both false negative and
false positive EWS [22, 23, 9]. Therefore, understanding the evolution of the
diffusion term $b(X_{t},t)\mathrm{d}W_{t}$ is crucial for reliable statements
on stability changes derived from data. Given that in real-world situations,
the assumption of time- and state-independent noise is hardly justifiable, a
more general theoretical framework advancing CSD to the case of time- and
state-dependent noise is called for.

In particular, a methodology is needed to extract from the observable a more
holistic picture of both the deterministic dynamics of the system and the
driving noise. The derivation of the variance and AC(1) in (3) hinges on the a
priori assumption of linear feedback. In applications, the system might
explore parts of the state space where non-linearities in the feedback are not
negligible anymore, putting the validity of Eq. (3) into question. In
contrast, the linear restoring rate $\lambda$ directly captures the desired
information of local stability and should therefore be considered the key
quantity to measure system stability and detect CSD. To obtain an estimation
$\widehat{\lambda}$, we perform a spatially local linear fit to the estimated
function $\widehat{a}(x)$ [24, 25] in some neighbourhood around
$\widehat{x}^{*}=\widehat{\mathbb{E}}[X]$ (see Methods and Supplementary
Material (SM) 1). A similar approach has recently been proposed in [26]. We
carry the concept to multiple dimensions and include an estimation of the
diffusion matrix $BB^{\top}(x,t)$ to supplement the standard CSD indicators
and to avoid false positives and false negatives caused by changes in the
driving noise. In particular, we discuss situations where the conventional CSD
indicators give an ambiguous or misleading picture. We show how the method
proposed herein conclusively resolves these ambiguities.

##  II Methods

It can be shown that the drift and diffusion coefficients $A$ and $B$ have the
following representation in terms of the increments
$\Delta\mathbf{X}_{t}:=\mathbf{X}_{t+\Delta t}-\mathbf{X}_{t}$ of the process
$\mathbf{X}$ [27, 28, 29, 24, 30, 31, 32, 33]:

$\displaystyle A(\mathbf{x},t)$ $\displaystyle=\lim_{\Delta t\rightarrow
0}\frac{1}{\Delta
t}\mathbb{E}\left[\Delta\mathbf{X}_{t}|\mathbf{X}_{t}=\mathbf{x}\right],$ (6)
$\displaystyle BB^{\top}(\mathbf{x},t)$ $\displaystyle=\lim_{\Delta
t\rightarrow 0}\frac{1}{\Delta
t}\mathbb{E}\left[\Delta\mathbf{X}_{t}\Delta\mathbf{X}_{t}^{\top}|\mathbf{X}_{t}=\mathbf{x}\right].$
(7)

If the stochastic differential equation (SDE) (4) exhibits effective time
independence, i.e. $A(\mathbf{x},t)\equiv A(\mathbf{x})$ and
$B(\mathbf{x},t)\equiv B(\mathbf{x})$ in some observation time span and if the
sample path of $\mathbf{X}$ is available at sufficiently small time steps
$\Delta t>0$, one may estimate $A(\mathbf{x})$ and $BB^{\top}(\mathbf{x})$ by
replacing the above ensemble average by the mean of the observed increments.
The law of large numbers yields consistent estimators that converge to the
true $A$ and $BB^{\top}$, omitting here a small bias stemming from the non-
zero $\Delta t$.

We generalise the definition of local system stability $\lambda$ in Eq. (5)
from the one-dimensional setting to the multi-dimensional setting of Eq. (1).
Consider the Jacobian matrix of $A$ at an equilibrium point
$\mathbf{x}^{*}\in\mathbb{R}^{n}$ with $A(\mathbf{x}^{*},t)=\mathbf{0}$. Such
an equilibrium point is stable if and only if all eigenvalues
$(-\lambda_{k}+i\omega_{k})_{k=1,\dots,n}$ of the Jacobian matrix
$\mathrm{D}A(\mathbf{x}^{*},t)$ exhibit negative feedback $-\lambda_{k}<0$.
Accordingly, we regard the set of $(\lambda_{k})_{k=1,\dots,n}$ as a measure
of system stability.

We obtain an estimation of $A$ and $BB^{\top}$ from time series data using the
respective estimators given in [34, 25]. A multivariate ordinary least squares
regression is performed to extract an estimate of the matrix $\mathrm{D}A$
(see SM 1 for details). The real parts of the corresponding eigenvalues are
then assessed in their time evolution. In a windowed time series analysis,
negative trends in any of the $\lambda_{k}$ would indicate a destabilisation
of the equilibrium state, which may point to an upcoming abrupt transition.
The windows must be short enough to justify the required time independence of
the dynamics within each individual window and yet comprise a sufficient
amount of data. The windowed estimation of $A(\mathbf{x},t)$ and
$BB^{\top}(\mathbf{x},t)$ then reveals potential temporal changes in the
system’s stability.

We will show that the estimation of the diffusion matrix $BB^{\top}$ may in
certain situations explain peculiar behaviour in the conventional CSD
indicators, such as simultaneously decreasing variance and increasing AC1. If
instead of local stability, here measured by $\lambda$, one wishes to examine
mean exit times from stable equilibria [35], the estimated diffusion matrix
$BB^{\top}$ has additional important implications [36].

##  III Results

First, we will examine the merits of the proposed CSD indicator
$\widehat{\lambda}$ on the estimated drift coefficient in a conceptual, one-
dimensional system subjected to time-dependent noise. Second, we will adopt a
two-dimensional predator-prey model with state-dependent noise from the
literature and perform a multi-dimensional analysis.

###  III.1 Fold bifurcation with time-dependent noise

The prototypical structure employed for conceptualising abrupt transitions in
many natural systems is the fold bifurcation [2, 37]. Consider therefore the
SDE defined by

$\displaystyle\mathrm{d}X_{t}$
$\displaystyle=\left(-X_{t}^{2}+\alpha(t)\right)\\!\mathrm{d}t+\sigma(t)\mathrm{d}W_{t},$
(8)

where $\alpha$ is the bifurcation parameter and $\sigma$ the noise strength.
For positive $\alpha>0$, there exists a stable equilibrium at
$x^{*}(t)=\sqrt{\alpha}$ which vanishes at the critical threshold
$\alpha_{\text{crit}}=0$ (see Fig. 1a). We simultaneously ramp down the noise
strength $\sigma(t)$. Such an evolution should be understood as a change in
the nature of the omitted fast dynamics, which cannot be ruled out in many
applications [38, 39].

Figure 1: Application of CSD indicators for synthetic data generated by the
model (8) in the main text. (a) Sample paths for zero noise (red) and with
noise (black). The noise strength $\sigma$ is ramped linearly from $0.2$ to
$0.06$ over the integration time span. (b), (c) Conventional CSD indicators of
variance and lag-1 autocorrelations (black) are calculated on detrended
windows and plotted along with the theoretical values (red) obtained from the
time-local Ornstein–Uhlenbeck linearisation. The shaded bands represent the
68% confidence intervals on $N=1000$ samples. (d) Estimator
$\widehat{\lambda}$ as proposed in this work, along with the true value
$\lambda=-\partial_{x}a(x^{*}(t),t)$ (red). All estimations were performed on
running windows of length $10^{2}$ consisting of $10^{3}$ data points,
considering the sample time-step is $\Delta t=10^{-1}$. A traditional analysis
using the variance and AC(1) would lead to a missed alarm, given their
opposing trends; in contrast, the CSD methodology proposed here clearly
detects the forthcoming bifurcation and correctly attributes the negative
variance trend to the decreasing amplitude of the driving noise $\sigma(t)$.

The temporal evolution of the variance and AC(1) can be approximated by Eq.
(3) after linearising the system around the time-dependent equilibrium
$x^{*}(t)$ (red lines in Fig. 1b and c):

$\mathrm{Var}\left[X_{t}\right]\approx\frac{\sigma(t)^{2}}{4\sqrt{\alpha(t)}},\quad\mathrm{AC}_{X}(1)\approx\exp\left(-2\sqrt{\alpha(t)}\Delta
t\right).$ (9)

The time-dependent noise thus induces a deceiving downward trend in the
estimated variance of the system (Fig. 1b) alongside increasing AC(1) (Fig.
1c). The conflicting indications given by variance and AC(1) would mislead the
observer to conclude that no significant EWS is present. In contrast, the
estimation of the linearised feedback $\lambda$ (Fig. 1d) clearly indicates a
weakening of local system stability and thus the presence of CSD.

The apparently inconsistent results of the conventional CSD indicators can be
understood and reconciled by examining the structure of the drift and
diffusion coefficients $a(x,t)=x^{2}-\alpha(t)$ and $b(x,t)=\sigma(t)$ and
their evolution in time. The true quantities for the functions $a(x,t)$ and
$b(x,t)$ at different times $t$ are plotted in Fig. S1 in SM 1 along with the
estimations obtained during the procedure outlined in the methods section.
After disclosing the time dependency of the diffusion coefficient in Fig. S1c,
the diverging trends in variance and AC(1) can be correctly interpreted in a
CSD assessment. The decrease in variance can be attributed to the decreasing
noise strength, and the approach of a bifurcation can be confirmed. In an
analogous setting featuring no bifurcation but an increasing noise strength,
the incurred increase in variance can be attributed correctly and the false
alarm that a conventional CSD analysis would raise can be avoided.

The statistical quality of the estimator $\widehat{\lambda}$ with respect to
its distribution width is similar to that of the conventional indicators and
sufficiently good to ensure a high likelihood of a statistically significant
trend. This can be argued by checking that the confidence intervals at the
beginning and the end of the estimator time series do not overlap (see also SM
2).

###  III.2 Predator-prey model with state-dependent noise

Following the work of Bengfort et al. [40], we examine a predator-prey model
for oceanic plankton populations (see SM 3 for details). Since Bengfort et al.
consider this model under the assumption of no external disturbances in the
form of noise, we adopt a noise model from a related study [41]. Because
environmental variability usually does not influence the population sizes
directly, but rather their growth rates, a multiplicative noise term is often
employed [42, 43, 44, 45]. This leads us to investigate the following system
of SDEs

$\displaystyle\mathrm{d}P_{t}$
$\displaystyle=\xi^{-1}\left(rP_{t}\left(1-\frac{P_{t}}{K(\text{turb})}\right)-\frac{aP_{t}^{2}}{h(\text{turb})^{2}+P_{t}^{2}}Z_{t}\right)\mathrm{d}t$
(10) $\displaystyle\quad+\xi^{-1/2}\sigma_{P}P_{t}\mathrm{d}W_{t}^{P},$ (11)
$\displaystyle\mathrm{d}Z_{t}$
$\displaystyle=\left(\frac{aP_{t}^{2}}{h(\text{turb})^{2}+P_{t}^{2}}Z_{t}-mZ_{t}^{2}\right)\mathrm{d}t+\sigma_{Z}Z_{t}\mathrm{d}W_{t}^{Z},$
(12)

under the external forcing of ocean turbulence. Due to the quadratic mortality
term $mZ^{2}$ of the predator population $Z$, this system can exhibit multiple
stable equilibria and bifurcations as indicated in Fig. 2a. The two white
noise terms are assumed to be independent and their strengths $\sigma_{P}$ and
$\sigma_{Z}$ are chosen such that noise-induced tipping only occurs in close
proximity to the bifurcation point.

Here, we examine the performance of the conventional and the newly proposed
CSD indicators as the system approaches this bifurcation. Fig. 2a shows sample
paths for the predator and prey populations along with the stable and unstable
equilibria of the prey population $P$ as implied by the parameter value turb
at time $t$.

Figure 2: Application of the CSD indicators on time-series data obtained from
the predator-prey model of [40]. (a) Sample paths of the prey population $P$
(black) and the predator population $Z$ (grey). The stable and unstable
equilibria of $P$ in their dependence on $\text{turb}(t)$ are plotted in red.
(b) and (c) show the means of the conventional CSD indicators variance and
AC(1) over $N=1000$ samples. (d) Real parts $-\lambda_{1,2}$ corresponding to
the estimated eigenvalues of the local Jacobian matrix. The eigenvalues have
been assigned (by colour) to the two populations, as the corresponding
eigenspace basis aligns very well. All estimations were performed on running
windows of length $100$, meaning $5000$ data points at sampling rate $\Delta
t=2\cdot 10^{-2}$.

The most common approach to the assessment of CSD in a multi-dimensional
system such as this one is to first reduce the system to one dimension [8, 46,
47]. The centre manifold theorem states that in close proximity to a critical
bifurcation, the direction of lowest stability will be the one to experience
further destabilisation. For this reason, a principle component analysis is
often performed to determine a linear combination of system components that
exhibits the largest variance or AC(1) and can thus be suspected to be of the
lowest stability [16, 17, 21]. However, as can be seen in the example at hand,
away from the immediate proximity of the bifurcation point, the destabilising
direction need not be the direction of lowest stability. Here, the identified
direction of stability loss would be closely aligned with the predator
population $Z$, as it operates on a slower time scale. This is problematic, as
this dimension is relatively impervious to changes in the control parameter
turb and will not exhibit CSD (grey curves in Fig. 2b and c).

To circumvent this issue, one should therefore perform a comprehensive
stability analysis on the multi-dimensional time series.

Figure 3: Prey population taken as a one-dimensional system (a) and
corresponding drift and diffusion coefficients $a(P,t)$ (b) and $b(P,t)$ (c)
for different time slices. In (b) and (c) the functions are plotted in their
$P$-dependence, while the time $t$ is represented by the respective colour of
the plot, as indicated in (a). The dashed lines in (b) represent the best
linear fits performed on the estimated $a(P,t)$. Annotated is the value of the
estimator $\hat{\lambda}$ on each window of data, i.e., the negative of the
slope of the respective linear fit. Similarly, linear fits in (c) are shown to
illustrate the apparent state-dependence, in addition to a linear fit of the
entire data (grey dashed).

This is achieved by examining both eigenvalues of the local equilibrium
dynamics as shown in Fig. 2d. The real part of the larger eigenvalue can be
seen to substantially decrease. This provides evidence for a destabilisation
along the more stable direction in the eigenspace. Note that the conventional
approach of focusing on the direction of largest eigenvalue would miss this
destabilisation.

Investigating again the conventional variance and AC(1) CSD indicators in Fig.
2b and c, AC(1) seems to indicate a destabilisation along the dimension of the
prey population $P$ but the trend of the observed variance seems to indicate
the opposite. This contradiction can be resolved by performing an analogous
analysis of the drift and diffusion coefficients for the one-dimensional time
series of $P$ as in the first example above (see Fig. 3). This reduction in
dimension can now be motivated by the fact that the vector in eigenspace
corresponding to the weakening eigenvalue in 2d lies predominantly in the
direction of $P$. The slopes of the estimated drift coefficient shown in Fig.
3b decrease as the system moves towards the bifurcation, agreeing with the
estimations in Fig. 2d. A clear state-dependence can be identified in the
estimated diffusion coefficient (Fig. 3c). Together with the observation of a
diminishing mean state in the prey population, it can be concluded that the
decrease in variance in 2b was due to a reduced noise amplitude and can thus
be reconciled with the increase in AC(1).

##  IV Discussion

Variance and AC(1) are often used in combination to assess whether or not a
system is approaching a critical transition. In general, a positive result is
considered robust when both indicators show a significant positive trend. We
have shown here that in the presence of time- or state-dependent noise
amplitudes, the variance of the system may actually decrease in the advent of
a bifurcation. If a monitored system shows a decreasing trend in variance
alongside an increasing trend in AC(1), this would typically not be considered
a robust EWS, leading to a missed alarm. An increase in noise strength over
time could, on the other hand, lead to a false alarm in form of an increasing
variance in systems with no underlying bifurcation. We have also shown that
common methods in dimension reduction can lead to missed alarms, as the
destabilising system component may not be the least stable to begin with.

To overcome these problems, we have proposed a method based on deriving a
Langevin equation from the observed dynamics. Our approach allows us to
separate the effects of possible CSD dynamics contained in the drift
coefficient from changes in the noise represented by the diffusion
coefficient. It also allows for a more holistic investigation of multi-
dimensional systems, without further mechanistic simplifications (see SM 4 for
a second example to this point, which shows in particular that the proposed
method works for periodic multi-dimensional systems that are problematic for
the conventional CSD indicators). We showed that our approach avoids the
pitfalls that a conventional CSD analysis suffers from in these examples.

We have shown that in the presented one-dimensional application, the
statistical quality of the estimator $\widehat{\lambda}$ is of the same order
as that of the estimators for variance and AC(1) (see also SM 2 for further
discussion). However, one important caveat bears mentioning: While for the
estimators of variance and AC(1), the length of the time series is the only
determining factor of convergence, the estimators for the drift and diffusion
coefficients also require small sampling time steps $1\gg\lambda\Delta t>0$ in
order for their bias to be small. In general, the estimator
$\widehat{\lambda}$ proposed here will still contain information about CSD
even in settings of large sample time steps $\Delta t$, but the signal-to-
noise ratio may prohibit its employment as a CSD indicator. Areas of
application where systems are potentially susceptible to tipping and where
high-frequency data may be available for analysis could be electricity grids
[48, 49, 50], financial markets [51, 52, 53], atmospheric circulation systems
such as monsoons [54, 55], ecosystems and vegetation systems such as the
Amazon rainforest, [56, 57, 58, 59], ocean circulation systems [3], or ice
sheets [60, 61].

Our method should be understood as a more general, reliable, and circumspect
indicator of CSD compared to the widely used variance and AC(1). Our approach
is appropriate in settings of generally time- and state-dependent driving
noise, where the combined conventional indicators fail. Moreover, the ability
to examine time series in their multi-dimensional complexity constitutes a
considerable improvement in the comprehension of the system compared to one-
dimensional summary statistics.

## Data Availability

Supplementary Material is available for download at the online version of this
manuscript. The implementation of the estimators introduced in this work is
available in the GitHub repository KramersMoyalEWS. Also included is the code
employed to generate all figures in the main text and the Supplementary
Material.

###### Acknowledgements.

This work has received funding from the Volkswagen Stiftung, the European
Union’s Horizon 2020 research and innovation programme under grant agreement
No. 820970 and under the Marie Sklodowska-Curie grant agreement No. 956170, as
well as from the Federal Ministry of Education and Research under grant No.
01LS2001A. This is TiPES contribution #X.

## Supplementary Material

### 1\. Details on the estimator $\widehat{\lambda}$

Here, we give a detailed description of the local stability measure
$\widehat{\lambda}$ in $n$-dimensional systems. For $n=1$, it is the negative
of the slope of the estimated drift coefficient $a$ around the equilibrium,
and therefore always a real number. For $n>1$, the estimation procedure
returns $n$ eigenvalues of the local Jacobian matrix of the estimated drift
coefficient in their algebraic multiplicity. The eigenvalues may be complex,
and thus the value of interest investigated in the main text is the negative
of the real part of the respective eigenvalues.

Figure S1: Estimated drift and diffusion coefficients $a(X,t)$ and $b(X,t)$
[25] for the one-dimensional fold bifurcation example in the main text. The
functions are plotted in their $X$-dependence, while the time $t$ is
represented by the respective colour of the plot. The dashed lines show the
theoretical quantities.

Prior to any analysis, the windowed time series data of each of the $n$
dimensions is linearly detrended. For the assessment of the conventional CSD
indicators, the mean of the data is removed, as they rely purely on the
fluctuations around the equilibrium state. In contrast, drift and diffusion
are assessed without subtraction of the mean to retain information about the
corresponding state dependence. In order to obtain numerical stability, the
$n$ time series are normalised to a standard deviation of 1, with no
implication on the subsequent estimations. For each window, the estimation of
the function $A(\mathbf{x})$ is returned as an array of values

$(\widehat{A}(\mathbf{x}_{i}))_{i=1,\dots,M^{n}},$ (13)

where $M$ is the number of evenly spaced target bins in each dimension. The
$M^{n}$ target bins in $\mathbb{R}^{n}$, therefore, form a grid on the
hypercube spanned by the state space explored by the time series. For each of
these bins, an estimation $\widehat{A}(\mathbf{x}_{i})$ is calculated using an
Epanechnikov kernel

$K(\mathbf{x})=\frac{3}{4h}\left(1-\frac{||\mathbf{x}||^{2}}{h^{2}}\right),~{}\mathrm{with}~{}\mathrm{support}~{}||\mathbf{x}||<h,$
(14)

with a kernel bandwidth $h$ of $14n/M$. Since the estimator
$\widehat{A}(\mathbf{x}_{i})$ will converge to some (biased) value as the
number of samples $\mathbf{X}_{k\Delta t}$ in the bin $\mathbf{x}_{i}$ tends
to infinity, it is clear that those bins with many samples converge fastest.
In our setting with equilibrium dynamics around one stable equilibrium
$\mathbf{x}^{*}$, this means that the estimations for bins closest to
$\mathbf{x}^{*}$ converge fastest, and the quality deteriorates for outer
bins. For this reason and in order to curtail the effects of a non-linear
drift term, we opt to only carry $50\%$ of bins centred around the bin
containing $\hat{\mathbf{x}}^{*}=\widehat{\mathbb{E}}[\mathbf{X}]$ to the
subsequent analysis. This is to say, we select a hypercube with side lengths
50% as large as the original hypercube. Thus, we are confronted with fixing
three free parameters a priori: The number of total bins $M$ in each
dimension, the percentage $m$ of bins to carry on either side of the estimated
equilibrium $\widehat{\mathbf{x}}^{*}$, and the kernel bandwidth. In this
study, we chose $M=50$ and $m=50\%$, meaning that for $n=1$, we have $25$
relevant bins for further analysis. The bandwidth is chosen as a function of
$M$ and $n$, as described above. However, the performance of the estimator
$\widehat{\lambda}$ is not very sensitive to small changes in these
parameters.

To obtain an estimation of the local Jacobian matrix around the equilibrium
point $\widehat{\mathbf{x}}^{*}$, we perform a multivariate ordinary linear
regression between $(\widehat{A}(\mathbf{x}_{i}))$ and $(\mathbf{x}_{i})$ over
$i$, including an intercept in the design matrix. The algebraic eigenvalues of
the resulting matrix are computed numerically. For $n=1$, this procedure is
equivalent to finding a best linear fit $(c-\widehat{\lambda}x_{i})$ to
$(\widehat{a}(x_{i}))$ over $i$.

### 2\. Assessing the statistical quality of $\hat{\lambda}$

The two applications presented in the main text demonstrate that CSD manifests
itself in a substantial negative trend of the estimator $\hat{\lambda}$ when
enough data are available. In this section, we aim to make this statement
concrete and to compare the indicator’s performance to that of variance and
AC(1). The assessment is based on the width of the three indicators’ numerical
distribution after application to synthetic data in one dimension. The top row
of Fig. S2 shows the distributions of the estimators that arise from the
application of the indicators to 1000 synthetic time series generated by
numerically integrating a time-homogeneous OU process:

$\mathrm{d}X_{t}=-\lambda X_{t}\mathrm{d}t+\mathrm{d}W_{t},$ (15)

Figure S2: Distributions of (a) variance, (b) AC(1), and (c) the estimator
$\hat{\lambda}$ on $1000$ samples of Ornstein–Uhlenbeck time series of length
$T=100$ each (see equation (15)). The sampling time step was set at $\Delta
t=0.1$. This estimation setting is equivalent to that chosen for the synthetic
fold bifurcation example in the main text. The blue histogram shows the
estimator distribution for the OU-parameter $\lambda=1$ and the violet for
$\lambda=0.1$. If the $97.5$\- and $2.5$-percentile of the two distributions
do not overlap, it is highly likely that the corresponding CSD will be
detected by the indicators. The separation in distribution is investigated for
several choices of $T$ and $\Delta t$ in the bottom row of the figure. Green
indicates no significant overlap of the estimator distributions, i.e., a
successful CSD estimation; red the contrary.

To analyse the behaviour of the estimators in a generic CSD scenario, i.e., a
temporal reduction of the restoring rate, we plot their distributions for
$\lambda=1$ and $\lambda=0.1$. If the distributions are sufficiently distinct,
the indicators may correctly detect a given reduction of the restoring rate
with a high likelihood. For different choices of window lengths and time steps
$\Delta t$, we check numerically whether this condition is satisfied for each
estimator (bottom row of Fig. S2). Being more sensitive at low data
availability, the estimators for variance and AC(1) perform better than that
for $\hat{\lambda}$. Above a window length of $T=100$, this difference is
negligible, judging by the proposed metric. The difference may also be less
pronounced when performing the same test on time series generated by models
with non-linear drift or jump-noise, where the state-locality of our method
can alleviate non-linear effects on the far ends of the state space.
Therefore, in a large range of applications, the estimator $\hat{\lambda}$
offers a statistically equally performant method of assessing CSD with the
additional advantage of robustness with respect to time- and state-dependent
noise, substantial advantages in higher-dimensional settings, as well as
settings featuring non-linear drifts and jumps in the noise.

Table SI: Parameter values used in the simulation of plankton populations following the model in equations (20). parameter | value | description
---|---|---
$r$ | $1$ | growth rate factor of prey $P$
$a$ | $1/9$ | rate of the predator consuming the prey
$m$ | $0.0525$ | mortality rate of the predator
$\xi$ | $0.7$ | time scale separation between prey and predator evolutions
$h_{0}$ | $1/16$ | factor influencing maximal consumption at zero turbulence
$c_{h}$ | $0.88$ | linear relationship between turbulence and $h$
$K_{0}$ | $0.7$ | carrying capacity at zero turbulence
$c_{K}$ | $0.66$ | linear relationship between turbulence and $K$
$\sigma_{P}$ | $0.037$ | strength of noise coupling to $P$
$\sigma_{Z}$ | $0.01$ | strength of noise coupling to $Z$

### 3\. Details on the predator-prey model

The specific model introduced in the main text is a modification of the
Truscott–Brindley model for ocean plankton populations originally introduced
in [62]. Bengfort et al. [40] generalised the model by introducing the
environmental parameter of fluid turbulence to the system and allowing higher
powers in the mortality term of the predator population. The full system
equations are given by

$\displaystyle\xi\dot{p}(t)$
$\displaystyle=rp(t)\left(1-\frac{p(t)}{K(\text{turb})}\right)-\frac{ap(t)^{2}}{h(\text{turb})^{2}+p(t)^{2}}z(t)$
(16) $\displaystyle\dot{z}(t)$
$\displaystyle=\frac{ap(t)^{2}}{h(\text{turb})^{2}+p(t)^{2}}z(t)-mz(t)^{2}$
(17) $\displaystyle K(\text{turb})$
$\displaystyle=K_{0}+c_{K}\cdot\text{turb}$ (18) $\displaystyle
h(\text{turb})$ $\displaystyle=\frac{h_{0}}{1+c_{h}\cdot\text{turb}}$ (19)

This system has been non-dimensionalised in order to reduce the number of
parameters. However, to retrieve realistic values of population sizes in units
of density, $p$ and $z$ merely need to be multiplied with constants $p_{0}$
and $z_{0}$. The first term on the right-hand side of the prey population’s
evolution $\dot{p}$ is the population growth rate as determined by the
relationship between the current population size and the carrying capacity
$K$. Below that capacity, the population grows and vice versa. The second term
is the mortality rate of the prey population, which is simultaneously the
growth rate of the predator population since it is assumed that all death in
$p$ and growth in $z$ occurs through consumption of the former by the latter.
The second term in the evolution of $z$ in the second equation is the
quadratic mortality term alluded to in the main text. This ultimately
facilitates multiple stable states as opposed to the same model with a linear
mortality term. The turbulence $\text{turb}\in[0,1]$ describes the normalised
strength of spatial mixing in the ocean modelled by circular eddies. All
parameter values but those for $\xi$, $c_{K}$ and $c_{h}$ are adopted directly
from [40] and can be found in Table SI along with a short description of their
interpretation. The parameters $c_{K}$ and $c_{h}$ were increased by a factor
of $2.2$ each for the purposes of this study to facilitate a bigger range of
stable prey populations in the large population regime. The fundamental nature
of the model remains unaltered by this change. Lastly, as described in the
main text, we introduced multiplicative noise terms commonly used in the
relevant literature [42, 43, 44, 45] to model environmental impacts on the
growth and mortality rates of the two populations. This leads us to the
complete set of model equations:

$\displaystyle\mathrm{d}P_{t}$
$\displaystyle=\xi^{-1}\\!\left(\\!rP_{t}\\!\left(1\\!-\\!\frac{P_{t}}{K(\text{turb})}\right)\\!-\\!\frac{aP_{t}^{2}}{h(\text{turb})^{2}\\!+\\!P_{t}^{2}}Z_{t}\\!\right)\\!\mathrm{d}t$
$\displaystyle\quad+\xi^{-1/2}\sigma_{P}P_{t}\mathrm{d}W_{t}^{P}.$
$\displaystyle\mathrm{d}Z_{t}$
$\displaystyle=\left(\frac{aP_{t}^{2}}{h(\text{turb})^{2}\\!+\\!P_{t}^{2}}Z_{t}-mZ_{t}^{2}\\!\right)\\!\mathrm{d}t+\sigma_{Z}Z_{t}\mathrm{d}W_{t}^{Z}.$
$\displaystyle K(\text{turb})$ $\displaystyle=K_{0}+c_{K}\cdot\text{turb},$
(20) $\displaystyle h(\text{turb})$
$\displaystyle=\frac{h_{0}}{1+c_{h}\cdot\text{turb}}$
$\displaystyle\text{turb}(t)$ $\displaystyle=1-\frac{7}{10}\frac{t}{T}.$

### 4\. Additional example of the multi-dimensional stability analysis

Figure S3: Employing the outlined method for multidimensional stability
analysis. (a) Time series data generated by a model undergoing a subcritical
Hopf bifurcation in normal form. (b) Estimates of the stability indicators
$\lambda_{k}$. The thick line represents the mean of the 200 sample
estimations. Since the Jacobian matrices at each point in time have complex
eigenvalues corresponding to the oscillatory dynamics, the eigenvalues are
conjugates and coincide in their real part. The theoretical value for the
negative of this real part is plotted in red.

In the example of the two-dimensional predator-prey model in the main text, it
was revealed by the analysis of the two-dimensional drift coefficient that the
dynamics could also be well-represented by two uncoupled one-dimensional SDEs
for the predator and prey population, respectively. As a result, the CSD
analysis was also comprehensive after a reduction to the prey dimension.
However, the dynamics of many systems cannot be reduced in such a way. This is
especially relevant for systems exhibiting pronounced oscillations. Using the
method outlined above, we therefore additionally assess the local stability of
a system undergoing a subcritical Hopf bifurcation in normal form.

$\displaystyle\mathrm{d}\mathbf{X}_{t}$
$\displaystyle=\begin{pmatrix}-\left(\mu(t)-\left(\mathbf{X}^{(1)}\right)^{2}-\left(\mathbf{X}^{(2)}\right)^{2}\right)\mathbf{X}^{(1)}-\omega\mathbf{X}^{(2)}\\\
-\left(\mu(t)-\left(\mathbf{X}^{(1)}\right)^{2}-\left(\mathbf{X}^{(2)}\right)^{2}\right)\mathbf{X}^{(2)}+\omega\mathbf{X}^{(1)}\end{pmatrix}\mathrm{d}t$
(21) $\displaystyle\quad+\varepsilon\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\mathrm{d}\mathbf{W}_{t},$ (22)

where $\omega=1$, $\varepsilon=0.01$, and $\mu(t)$ decreases linearly from $2$
to $0.1$ over the integration time of $T=1000$. For $\mu>0$, the origin is a
stable fixed point with eigenvalues $-\mu\pm i\omega$. Furthermore, there is
an unstable limit cycle with radius $\sqrt{\mu}$ and perturbations from the
origin decay in the form of spirals. At $\mu=0$, the radius of the unstable
limit cycle reaches zero, and the origin turns into an unstable fixed point.
The data was sampled at time steps $\Delta t=0.1$ and analysed in windows of
length $T=100$. The results are presented in Fig. S3. The real parts of both
eigenvalues are known to be $-\mu(t)$, and the estimations track this value
relatively closely. A destabilisation of the equilibrium can clearly be made
out. An additional insight gained via the CSD assessment through the Langevin
equation approach proposed here is that the local system exhibits oscillatory
dynamics, as identified by the complex eigenvalues of the Jacobian matrix.

## References

  * Lenton  _et al._ [2008] T. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf, and H. J. Schellnhuber, Tipping elements in the Earth’s climate system, Proceedings of the National Academy of Sciences of the United States of America 105, 1786 (2008).
  * Scheffer  _et al._ [2009] M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. van Nes, M. Rietkerk, and G. Sugihara, Early-warning signals for critical transitions, Nature 461, 53 (2009).
  * Boers [2021] N. Boers, Observation-based early-warning signals for a collapse of the Atlantic Meridional Overturning Circulation, Nature Climate Change 11, 680 (2021).
  * Boers  _et al._ [2022] N. Boers, M. Ghil, and T. F. Stocker, Theoretical and paleoclimatic evidence for abrupt transitions in the earth system, Environmental Research Letters 17, 093006 (2022).
  * Zwanzig [2001] R. Zwanzig,  _Nonequilibrium Statistical Mechanics_ (Oxford University Press, 2001).
  * Hasselmann [1976] K. Hasselmann, Stochastic climate models Part I. Theory, Tellus 28, 473 (1976).
  * Horsthemke and Lefever [1984] W. Horsthemke and R. Lefever,  _Noise-Induced Transitions_ (Springer Berlin, 1984).
  * Meisel  _et al._ [2015] C. Meisel, A. Klaus, C. Kuehn, and D. Plenz, Critical slowing down governs the transition to neuron spiking, PLOS Computational Biology 11, e1004097 (2015).
  * Boettner and Boers [2022] C. Boettner and N. Boers, Critical slowing down in dynamical systems driven by nonstationary correlated noise, Physical Review Research 4, 013230 (2022).
  * Morr  _et al._ [2022] A. Morr, D. Kreher, and N. Boers, Red noise in continuous-time stochastic modelling, arXiv (2022).
  * Ashwin  _et al._ [2012] P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox, Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 1166 (2012).
  * Dakos  _et al._ [2008] V. Dakos, M. Scheffer, E. H. van Nes, V. Brovkin, V. Petoukhov, and H. Held, Slowing down as an early warning signal for abrupt climate change, Proceedings of the National Academy of Sciences 105, 14308 (2008).
  * deMenocal  _et al._ [2000] P. deMenocal, J. Ortiz, T. Guilderson, J. Adkins, M. Sarnthein, L. Baker, and M. Yarusinsky, Abrupt onset and termination of the African Humid Period: rapid climate responses to gradual insolation forcing, Quaternary Science Reviews 19, 347 (2000).
  * Boers  _et al._ [2017] N. Boers, N. Marwan, H. M. J. Barbosa, and J. Kurths, A deforestation-induced tipping point for the South American monsoon system, Scientific Reports 7, 41489 (2017).
  * Kuehn [2011] C. Kuehn, A mathematical framework for critical transitions: Bifurcations, fast–slow systems and stochastic dynamics, Physica D: Nonlinear Phenomena 240, 1020 (2011).
  * Held and Kleinen [2004] H. Held and T. Kleinen, Detection of climate system bifurcations by degenerate fingerprinting, Geophysical Research Letters 31, 10.1029/2004GL020972 (2004).
  * Dakos [2018] V. Dakos, Identifying best-indicator species for abrupt transitions in multispecies communities, Ecological Indicators 94, 494 (2018).
  * Kuehn [2013] C. Kuehn, A mathematical framework for critical transitions: Normal forms, variance and applications, Journal of Nonlinear Science 23, 457 (2013).
  * Lenton  _et al._ [2012] T. M. Lenton, V. N. Livina, V. Dakos, E. H. van Nes, and M. Scheffer, Early warning of climate tipping points from critical slowing down: comparing methods to improve robustness, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 1185 (2012).
  * Uhlenbeck and Ornstein [1930] G. E. Uhlenbeck and L. S. Ornstein, On the theory of the brownian motion, Phys. Rev. 36, 823 (1930).
  * Weinans  _et al._ [2021] E. Weinans, R. Quax, E. H. van Nes, and I. A. v. d. Leemput, Evaluating the performance of multivariate indicators of resilience loss, Scientific Reports 11, 9148 (2021).
  * Dakos  _et al._ [2012] V. Dakos, E. H. van Nes, P. D’Odorico, and M. Scheffer, Robustness of variance and autocorrelation as indicators of critical slowing down, Ecology 93, 264 (2012).
  * Clarke  _et al._ [2023] J. J. Clarke, C. Huntingford, P. D. L. Ritchie, and P. M. Cox, Seeking more robust early warning signals for climate tipping points: the ratio of spectra method (ROSA), Environmental Research Letters 18, 035006 (2023).
  * Lamouroux and Lehnertz [2009] D. Lamouroux and K. Lehnertz, Kernel-based regression of drift and diffusion coefficients of stochastic processes, Physics Letters A 373, 3507 (2009).
  * Rydin Gorjão and Meirinhos [2019] L. Rydin Gorjão and F. Meirinhos, kramersmoyal: Kramers–moyal coefficients for stochastic processes, Journal of Open Source Software 4, 1693 (2019).
  * Wand  _et al._ [2023] T. Wand, M. Heßler, and O. Kamps, Memory effects, multiple time scales and local stability in langevin models of the s&p500 market correlation (2023), arXiv:2307.12744 [q-fin.ST] .
  * Risken and Frank [1996] H. Risken and T. Frank,  _The Fokker–Planck equation_, 2nd ed. (Springer-Verlag, Berlin, Heidelberg, 1996).
  * van Kampen [2007] N. G. van Kampen, Chapter VIII - The Fokker–Planck equation, in  _Stochastic Processes in Physics and Chemistry_, North-Holland Personal Library (Elsevier, 2007) 3rd ed., pp. 193–218.
  * Gardiner [2009] C. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Sciences_ , 4th ed. (Springer-Verlag Berlin Heidelberg, 2009).
  * Honisch and Friedrich [2011] C. Honisch and R. Friedrich, Estimation of Kramers-Moyal coefficients at low sampling rates, Physical Review E 83, 066701 (2011).
  * Tabar [2019] M. R. R. Tabar,  _Analysis and Data-Based Reconstruction of Complex Nonlinear Dynamical Systems_, 1st ed. (Springer International Publishing, 2019).
  * Hassanibesheli  _et al._ [2020] F. Hassanibesheli, N. Boers, and J. Kurths, Reconstructing complex system dynamics from time series: a method comparison, New Journal of Physics 22, 073053 (2020).
  * Riechers  _et al._ [2023] K. Riechers, L. Rydin Gorjão, F. Hassanibesheli, P. G. Lind, D. Witthaut, and N. Boers, Stable stadial and interstadial states of the last glacial’s climate identified in a combined stable water isotope and dust record from greenland, Earth System Dynamics 14, 593 (2023).
  * Rydin Gorjão  _et al._ [2019] L. Rydin Gorjão, J. Heysel, K. Lehnertz, and M. R. R. Tabar, Analysis and data-driven reconstruction of bivariate jump-diffusion processes, Physical Review E 100, 062127 (2019).
  * Arani  _et al._ [2021] B. M. S. Arani, S. R. Carpenter, L. Lahti, E. H. van Nes, and M. Scheffer, Exit time as a measure of ecological resilience, Science 372, eaay4895 (2021).
  * Kramers [1940] H. A. Kramers, Brownian motion in a field of force and the diffusion model of chemical reactions, Physica 7, 284 (1940).
  * Guttal and Jayaprakash [2007] V. Guttal and C. Jayaprakash, Impact of noise on bistable ecological systems, Ecological Modelling 201, 420 (2007).
  * Fischer and Schär [2009] E. M. Fischer and C. Schär, Future changes in daily summer temperature variability: driving processes and role for temperature extremes, Climate Dynamics 33, 917 (2009).
  * Shiller [1988] R. J. Shiller, Causes of changing financial market volatility (1988) pp. 1–28, https://www.kansascityfed.org/Jackson%20Hole/documents/3811/1988-S88SHILL.pdf.
  * Bengfort  _et al._ [2014] M. Bengfort, U. Feudel, F. M. Hilker, and H. Malchow, Plankton blooms and patchiness generated by heterogeneous physical environments, Ecological Complexity 20, 185 (2014).
  * Sieber  _et al._ [2007] M. Sieber, H. Malchow, and L. Schimansky-Geier, Constructive effects of environmental noise in an excitable prey–predator plankton system with infected prey, Ecological Complexity 4, 223 (2007).
  * May [1973] R. M. May, Stability in randomly fluctuating versus deterministic environments, The American Naturalist 107, 621 (1973).
  * Samanta [2011] G. P. Samanta, A stochastic two species competition model: nonequilibrium fluctuation and stability, International Journal of Stochastic Analysis 2011, 489386 (2011).
  * Das and Samanta [2021] A. Das and G. P. Samanta, Influence of environmental noises on a prey–predator species with predator-dependent carrying capacity in alpine meadow ecosystem, Mathematics and Computers in Simulation 190, 1294 (2021).
  * Belabbas  _et al._ [2021] M. Belabbas, A. Ouahab, and F. Souna, Rich dynamics in a stochastic predator-prey model with protection zone for the prey and multiplicative noise applied on both species, Nonlinear Dynamics 106, 2761 (2021).
  * Boers [2018] N. Boers, Early-warning signals for Dansgaard-Oeschger events in a high-resolution ice core record, Nature Communications 9, 2556 (2018).
  * Bury  _et al._ [2020] T. M. Bury, C. T. Bauch, and M. Anand, Detecting and distinguishing tipping points using spectral early warning signals, Journal of The Royal Society Interface 17, 20200482 (2020).
  * Hines  _et al._ [2011] P. Hines, E. Cotilla-Sanchez, and S. Blumsack, Topological models and critical slowing down: Two approaches to power system blackout risk analysis, in  _2011 44th Hawaii International Conference on System Sciences_ (2011) pp. 1–10.
  * Cotilla-Sanchez  _et al._ [2012] E. Cotilla-Sanchez, P. D. H. Hines, and C. M. Danforth, Predicting critical transitions from time series synchrophasor data, IEEE Transactions on Smart Grid 3, 1832 (2012).
  * Ren and Watts [2015] H. Ren and D. Watts, Early warning signals for critical transitions in power systems, Electric Power Systems Research 124, 173 (2015).
  * Guttal  _et al._ [2016] V. Guttal, S. Raghavendra, N. Goel, and Q. Hoarau, Lack of critical slowing down suggests that financial meltdowns are not critical transitions, yet rising variability could signal systemic risk, PLOS ONE 11, e0144198 (2016).
  * Gatfaoui  _et al._ [2017] H. Gatfaoui, I. Nagot, and P. de Peretti, 3 - are critical slowing down indicators useful to detect financial crises?, in  _Systemic Risk Tomography_, edited by M. Billio, L. Pelizzon, and R. Savona (Elsevier, 2017) pp. 73–93.
  * Diks  _et al._ [2019] C. Diks, C. Hommes, and J. Wang, Critical slowing down as an early warning signal for financial crises?, Empirical Economics 57, 1201 (2019).
  * Levermann  _et al._ [2009] A. Levermann, J. Schewe, V. Petoukhov, and H. Held, Basic mechanism for abrupt monsoon transitions, Proceedings of the National Academy of Sciences 106, 20572 (2009).
  * Schewe  _et al._ [2012] J. Schewe, A. Levermann, and H. Cheng, A critical humidity threshold for monsoon transitions, Climate of the Past 8, 535 (2012).
  * Boulton  _et al._ [2022] C. A. Boulton, T. M. Lenton, and N. Boers, Pronounced loss of Amazon rainforest resilience since the early 2000s, Nature Climate Change 12, 271 (2022).
  * Smith  _et al._ [2022] T. Smith, D. Traxl, and N. Boers, Empirical evidence for recent global shifts in vegetation resilience, Nature Climate Change 12, 477 (2022).
  * Forzieri  _et al._ [2022] G. Forzieri, V. Dakos, N. G. McDowell, A. Ramdane, and A. Cescatti, Emerging signals of declining forest resilience under climate change, Nature 608, 534 (2022).
  * Smith and Boers [2023] T. Smith and N. Boers, Global vegetation resilience linked to water availability and variability, Nature Communications 14, 498 (2023).
  * Levermann and Winkelmann [2016] A. Levermann and R. Winkelmann, A simple equation for the melt elevation feedback of ice sheets, The Cryosphere 10, 1799 (2016).
  * Boers and Rypdal [2021] N. Boers and M. Rypdal, Critical slowing down suggests that the western greenland ice sheet is close to a tipping point, Proceedings of the National Academy of Sciences 118, e2024192118 (2021).
  * Truscott and Brindley [1994] J. E. Truscott and J. Brindley, Ocean plankton populations as excitable media, Bulletin of Mathematical Biology 56, 981 (1994).


# Adversarial examples by perturbing high-level features
in intermediate decoder layers

Vojtěch Čermák1, Lukáš Adam1

###### Abstract

We propose a novel method for creating adversarial examples. Instead of
perturbing pixels, we use an encoder-decoder representation of the input image
and perturb intermediate layers in the decoder. This changes the high-level
features provided by the generative model. Therefore, our perturbation
possesses semantic meaning, such as a longer beak or green tints. We formulate
this task as an optimization problem by minimizing the Wasserstein distance
between the adversarial and initial images under a misclassification
constraint. We employ the projected gradient method with a simple inexact
projection. Due to the projection, all iterations are feasible, and our method
always generates adversarial images. We perform numerical experiments on the
MNIST and ImageNet datasets in both targeted and untargeted settings. We
demonstrate that our adversarial images are much less vulnerable to
steganographic defence techniques than pixel-based attacks. Moreover, we show
that our method modifies key features such as edges and that defence
techniques based on adversarial training are vulnerable to our attacks.

##  1 Introduction

In the past decade, the widespread application of deep neural networks raised
security concerns as it creates incentives for attackers to exploit any
potential weakness. For example, attackers could create traffic signs
invisible to autonomous cars or make malware filters ignore threats. Those
security concerns rose since (Szegedy et al. 2014) showed that deep neural
networks are vulnerable to small perturbations of inputs that are designed to
cause misclassification of a classifier. These adversarial examples are
indistinguishable from natural examples as the perturbation is too small to be
perceived by humans.

We extend the usual approach to constructing adversarial examples that focuses
on norm-bounded pixel modifications. Instead of perturbing the pixels
directly, we perturb features learned by an encoder-decoder model. There are
several ways to perturb the features from the decoder, ranging from high-level
features collected from initial decoder layers to much finer features
collected from decoder layers near the reconstructed image. When we perturb
features in the initial layers near the latent image representation, even
small perturbations may completely change the image meaning. On the other
hand, perturbations on fine features are very close to pixel perturbations and
carry little semantic information. As a compromise, we suggest perturbing
intermediate decoder layers.

SourceLatentIntermediateReconstructedEncoderDecoder1Decoder2intermediate
latent representationinserted perturbation

Our way of generating adversarial images provides several advantages:

  * •

The perturbations are interpretable. They modify high-level features such as
fur colour or beak length.

  * •

The perturbations often follow edges. This means that they are less detectable
by steganographic defence techniques (Johnson and Jajodia 1998).

  * •

The perturbations keep the structure of the unperturbed images, including
hidden structures unrelated to their semantic meaning.

  * •

The decoder ensures that the adversarial images have the correct pixel values;
there is no need to project to $[0,1]$.

The drawback of our approach is that we cannot perturb an arbitrary image but
only those representable by the decoder.

To obtain a formal optimization problem, we minimize the distance between the
original and adversarial images in the reconstructed space. To successfully
generate an adversarial image, we add a constraint requiring that the
perturbed image is misclassified. To handle this constraint numerically, we
propose to use the projected gradient method. We approximate the projection
operator by a fast method which always returns a feasible point. This brings
additional benefits to our method:

  * •

The method works with feasible points. Even if it does not converge, it
produces an adversarial image.

  * •

Because we minimize the distance between the original and adversarial images,
we do not need to specify their maximal possible distance as many methods do.

Numerical experiments present results on MNIST and ImageNet in both targeted
and untargeted settings. They support all advantages mentioned above. We use
both the $l_{2}$ and Wasserstein distances and show in which settings each
performs better. We select a simple steganography defence mechanism and show
that our attack is much less vulnerable to it than pixel-based perturbations.
We examine the difference between perturbing the latent and intermediate
layers and show that the intermediate layers indeed modify the high-level
features of the reconstructed image. Finally, we show how our algorithm
gradually incorporates the high-level features from the initial to the
adversarial image. Our codes are available online to promote
reproducibility.111https://anonymous.4open.science/r/latent-adv-examples-5C59

### Related work

The threat model based on small $l_{p}$ norm-bounded perturbations has been
the main focus of research. It was originally introduced in (Szegedy et al.
2014), where they used L-BFGS to minimize the $l_{2}$ distance between the
original and adversarial images. Since then, many new ways how to both attack
and defend against adversarial examples have been introduced. (Goodfellow,
Shlens, and Szegedy 2015) introduced Fast Gradient Sign Method (FGSM), which
is bounded by the $l_{\infty}$ metric. The authors used FGSM to generate new
adversarial examples and used them to augment the training set in the
adversarial training defence technique. A natural way to extend the FGSM
attack is to iterate the gradient step as it is done in the Basic Iterative
Method of (Kurakin, Goodfellow, and Bengio 2016) and Projected Gradient
Descend attack of (Madry et al. 2018). (Papernot et al. 2016) argued to use
the $l_{0}$ distance to model human perception and proposed a class of attacks
optimized under $l_{0}$ distance. (Carlini and Wagner 2017) designed strong
attack algorithms based on optimizing the $l_{0}$, $l_{2}$ and $l_{\infty}$
distances.

(Xiao et al. 2018) used GANs to generate noise for adversarial perturbation.
Other authors used generative models in defence against adversarial examples.
MagNet of (Meng and Chen 2017) used reconstruction error of variational
autoencoders to detect adversarial examples. A similar idea was used in
DefenceGAN of (Samangouei, Kabkab, and Chellappa 2018), where they cleaned
adversarial images by matching them to their representation in latent space of
the GAN trained on clean data.

We use similar optimization techniques as in decision-based attacks such as
Boundary attack of (Brendel, Rauber, and Bethge 2018) and HopSkipJump attack
of (Chen, Jordan, and Wainwright 2020). The optimization techniques always
keep intermediate results in the adversarial region and always outputs
misclassified examples.

Some works have already investigated adversarial examples outside of the
$l_{p}$ norm. (Wong, Schmidt, and Kolter 2019) used the Sinkhorn approximation
of the Wasserstein distance to create adversarial images with the same
geometrical structure as the original images. The Wasserstein distance is more
suitable than the standard $l_{p}$ distance metrics because it preserves
differences in high-level features. We increase this benefit by additionally
modifying the high-level features instead of pixels.

The Unrestricted Adversarial Examples of (Song et al. 2018), are adversarial
examples constructed from scratch using conditional generative models. Our
paper differs in several aspects. First, we use the intermediate instead of
the latent representation to perturb high-level features. Second, we use an
unconditional generator, which allows us to freely move in the intermediate
latent space and perform several operations such as the projection.

##  2 Proposed Formulation

Finding an adversarial image amounts to finding some image $x$ which is close
to a given $x_{0}$, and the neural network misclassifies it. For reasons
mentioned in the introduction, we do not work with the original images but
with their latent representations. Therefore, we need a decoder $D$ which maps
the latent representation $z$ to the original representation $x$. Since we
want to insert some perturbation $p$ to the intermediate layer in the decoder,
we split the decoder into two parts $D=D_{2}\circ D_{1}$. We call $D_{1}(z)$
the intermediate latent representation and $D_{2}(D_{1}(z))$ the reconstructed
image.

The mathematical formulation of finding an adversarial image close to $x_{0}$
reads:

$\displaystyle\operatorname*{minimize}_{p}$
$\displaystyle\operatorname{distance}(x,x_{0})$ (1) subject to $\displaystyle
x=D_{2}(D_{1}(z_{0})+p),$ $\displaystyle x_{0}=D_{2}(D_{1}(z_{0})),$
$\displaystyle g(x)\leq 0,$ $\displaystyle x\in[0,1]^{n}.$

Here, $z_{0}$ is the latent representation of the image $x_{0}$. We insert the
perturbation $p$ to the intermediate latent representation $D_{1}(z_{0})$ and
only then reconstruct the image by applying the second part of the decoder
$D_{2}$.

When $D_{1}$ is the identity, and $D_{2}$ is the decoder, we obtain the case
when perturbations appear in the latent space. Similarly, when $D_{1}$ is the
decoder, and $D_{2}$ is the identity, we recover the standard pixel
perturbations. Therefore, our formulation (1) generalizes both approaches.

### Objective

The objective function measures the distance between the adversarial $x$ and
the original $x_{0}$ image in the reconstructed space. We will use the $l_{2}$
norm and the Wasserstein distance. Having the distance in the objective is
advantageous because we do not need to specify the $\varepsilon$-neighborhood
when this distance is in the constraint.

### Constraints

Formulation (1) has multiple constraints. Constraint $x\in[0,1]^{n}$ says that
the pixels need to stay in this range. Since $x$ is an output of the decoder,
this constraint is always satisfied.

The other constraint $g(x)\leq 0$ is a misclassification constraint. We employ
the margin function

$m(x,k)=\underset{i,i\neq k}{\text{max }}F_{i}(x)-F_{k}(x),$ (2)

where $F=[F_{1}(x),...,F_{m}(x)]$ is a classifier with $m$ classes and $k$ is
a class index. We define this constraint for both targeted and untargeted
attacks:

$\displaystyle\text{targeted}:$ $\displaystyle g(x)=m(x,k),$ (3)
$\displaystyle\text{untargeted}:$ $\displaystyle
g(x)=-m(x,\operatorname{argmax}_{i=1,\dots,m}F_{i}(x)).$

For the targeted case, we need to specify the target class $k$, while for the
untargeted case, the class $k$ is the classifier prediction. For the latter
case, we need to multiply the margin by $-1$ as we require the misclassified
images to satisfy $g(x)\leq 0$.

We will later use that the original image always satisfies $g(x_{0})>0$
because otherwise the original image is already misclassified, and no
perturbation $(p=0)$ is the optimal solution of (1).

##  3 Proposed Solution Method

We propose to solve (1) by the projected gradient method (Nocedal and Wright
2006) with inexact projection. Our algorithm produces a feasible point at
every iteration. Therefore, it always generates a misclassified image.

### Proposed algorithm

Since we optimize with respect to $p$, we define the objective and constraint
by

$\displaystyle f(p)$
$\displaystyle=\operatorname{distance}(D_{2}(D_{1}(z_{0})+p),x_{0}),$ (4)
$\displaystyle\hat{g}(p)$ $\displaystyle=g(D_{2}(D_{1}(z_{0})+p)).$

We describe our procedure in Algorithm 1. First, we initialize $p$ by some
feasible $p_{\rm init}$ and then run the projected gradient method for
$\max_{\rm iter}$ iterations. Step 8 computes the optimization step by
minimizing the objective $f$ and step 9 uses the inexact projection (described
later) to project the suggested iteration $p_{\rm next}$ back to the feasible
set. As we will see later, the constraint $p_{\rm next}\neq p$ implies that
$p_{\rm next}$ was not feasible, and it was projected onto the boundary. In
such a case, steps 4-7 “bounce away” from the boundary. Since
$-\nabla\hat{g}(p)$ points inside the feasible set, step 5 finds some
$\beta>0$ such that $p-\beta\nabla\hat{g}(p)$ lies in the interior of the
feasible set.

Algorithm 1 for finding adversarial images by solving (1)

1:$p\leftarrow p_{\rm init}$

2:for $i\in\\{0,\dots,\max_{\rm iter}\\}$ do

3: if $i>0$ and $p_{\rm next}\neq p$ then

4: $\beta\leftarrow\beta_{i}$

5: while $\hat{g}(p-\beta\nabla\hat{g}(p))\geq 0$ do

6: $\beta\leftarrow\frac{\beta}{2}$

7: $p\leftarrow p-\beta\nabla\hat{g}(p)$

8: $p_{\rm next}\leftarrow p-\alpha_{i}\nabla f(p)$

9: $p\leftarrow\textsc{Project}(p_{\rm next},p,\delta)$

10:return $p$

### Inexact projection

Step 9 in Algorithm 1 uses the inexact projection. We summarize this
projection in Algorithm 2. Its main idea is to find a feasible point on the
line between $p$ and $p_{\rm next}$. This line is parameterized by
$(1-c)p+cp_{\rm next}$ for $c\in[0,1]$. Due to the construction of Algorithm
1, $p$ is always strictly feasible and therefore $\hat{g}(p)<0$. If
$\hat{g}(p_{\rm next})<0$, then $p_{\rm next}$ is feasible and we accept it in
step 4. In the opposite case, we have $\hat{g}(p)<0$ and $\hat{g}(p_{\rm
next})\geq 0$ and we can use the bisection method to find some $c\in(0,1)$ for
which the constraint is satisfied. The standard bisection method would return
a point with $\hat{g}((1-c)p+cp_{\rm next})=0$, however, since the formulation
(1) contains the inequality constraint, we require the constraint value to lie
only in some interval $[-\delta,0]$ instead of being zero.

Algorithm 2 Inexact projection onto the feasible set

1:procedure Project($p_{\rm next},p,\delta$)

2: assert $\hat{g}(p)<0$

3: if $\hat{g}(p_{\rm next})<0$ then

4: return $p_{\rm next}$

5: $a\leftarrow 0$, $b\leftarrow 1$, $c\leftarrow\frac{1}{2}(a+b)$

6: while $\hat{g}((1-c)p+cp_{\rm next})\notin[-\delta,0]$ do

7: if $\hat{g}((1-c)p+cp_{\rm next})>0$ then

8: $b\leftarrow c$

9: else

10: $a\leftarrow c$

11: $c\leftarrow\frac{1}{2}(a+b)$

12: return $(1-c)p+cp_{\rm next}$

Figure 1 depicts the projection in the intermediate latent space. The light-
grey region is feasible, and the white region is infeasible. The infeasible
region always contains $D_{1}(z_{0})$, while the feasible region always
contains $D_{1}(z_{0})+p$. If $D_{1}(z_{0})+p_{\rm next}$ lies in the feasible
region, the projection returns $p_{\rm next}$. In the opposite case, we look
for some point on the line between $D_{1}(z_{0})+p$ and $D_{1}(z_{0})+p_{\rm
next}$. The points which may be accepted are depicted by the thick solid line.
The length of this line is governed by the threshold $\delta$. The extremal
case $\delta=0$ always returns the point on the boundary, while
$\delta=\infty$ prolongs the line to $D_{1}(z_{0})+p$.

Figure 1: Inexact projection in the intermediate latent space. The points
which may be accepted are depicted by the thick solid line.Feasible
regionInfeasible region$D_{1}(z_{0})+p_{\rm
next}$$D_{1}(z_{0})+p$$D_{1}(z_{0})$ Accepted projection

### Analysis of the algorithm

The preceding text mentioned that Algorithm 1 generates a strictly feasible
point $p$, thus $g(D_{2}(D_{1}(z_{0})+p))<0$. We prove this in the next
theorem and add a speed of convergence. We postpone the proof to the appendix.

######  Theorem 1.

Assume that $g\circ D_{2}$ is a continuous function. Then Algorithm 2
generates a strictly feasible point $p$ in a finite number of iterations. If
$g\circ D_{2}$ is moreover Lipschitz continuous with modulus $L$, then
Algorithm 2 converges in at most $\log_{2}\frac{L\|p_{\rm next}-p\|}{\delta}$
iterations.

The previous theorem implies that all iterations $p$ produced by Algorithm 1
are strictly feasible.

##  4 Numerical Experiments

This section shows numerical experiments on the image datasets MNIST (LeCun,
Cortes, and Burges 2010) and ImageNet (Deng et al. 2009).

### Used architectures

As the classifiers to be fooled in our MNIST experiments, we train a non-
robust classifier with architecture based on VGG blocks (Simonyan and
Zisserman 2015) and use the publicly available robust model (Madry et al.
2018). We generate new digits using an unconditional ALI generator (Donahue,
Krähenbühl, and Darrell 2016; Dumoulin et al. 2016). For ImageNet we use
EfficientNet B0 (Tan and Le 2019) as a classifier and BigBiGAN (Donahue and
Simonyan 2019) as an encoder-decoder model.

### Numerical setting

As an objective we use the standard $l_{2}$ distance and the Sinkhorn
approximation (Cuturi 2013) to the Wasserstein distance. This approximation
adds a weighted Kullback-Leibler divergence to solve

$\displaystyle\operatorname*{minimize}_{W}$ $\displaystyle\langle
C,W\rangle+\lambda\langle W,\log W\rangle$ (5) subject to $\displaystyle
W1=x_{0},W^{\top}1=x,W\geq 0.$

The optimal value of this problem is the Sinkhorn distance between images
$x_{0}$ and $x$. The matrix $C$ specifies the distance between pixels. Since
$x_{0}$ and $x$ are required to be probability distributions, we normalize the
pixels from decoder output to sum to one. We do not need to impose the non-
negativity constraint because the decoder outputs positive pixel intensities.
We use the implementation from the GeomLoss package (Feydy et al. 2019). Since
the number of elements $W$ from (5) equals the number of pixels squared, using
the Sinkhorn distance was infeasible for ImageNet, where we use only the
$l_{2}$ distance.

For MNIST, we randomly generate the latent images $z_{0}$ and then use the
decoder for reconstructed images $x_{0}$. We required that the classifier
predicts the digit into the correct class with the probability of at least
0.99. For ImageNet, this technique produces images of lower quality, and we
fed the encoder with real images and used their encoded representation.

We run all experiments for 1000 iterations. Algorithm 1 requires stepsizes
$\alpha$ and $\beta$, while Algorithm 2 requires the threshold $\delta$. We
found that the stepsize $\alpha$ is not crucial because even if it is large,
the projection will reduce it. We therefore selected $\alpha=1$. For the
second stepsize $\beta$, we selected a simple annealing scheme with
exponential decay. In both cases, we used the normed gradient. For the
threshold we selected $\delta=1$. Since the margin function (2) is bounded by
$1$, Figure 1 then implies that the acceptable projection increases all the
way to $D_{1}(z_{0})+p$. Since at least one iteration of the projection is
performed, the distance to the boundary at least halves. We found this to be a
good compromise between speed and approximative quality. Theorem 1 then
implies that Algorithm 1 converges in at most $\log_{2}L$ iterations.

We run our experiments on an Nvidia GPU with 6GB of VRAM.

### Numerical results for MNIST

Network | Attack |  $l_{2}$ distance | Wasserstein distance
---|---|---|---
|  | CW attack |  $l_{2}$ attack | Wasserstein | CW attack |  $l_{2}$ attack | Wasserstein
Non-robust | Untargeted | 7.85 | 15.28 | 20.87 | 1.77 | 1.49 | 0.08
| Targeted | 10.07 | 16.78 | 24.01 | 2.76 | 1.83 | 0.12
Robust | Untargeted | 17.36 | 22.78 | 24.68 | 0.78 | 2.48 | 0.19
| Targeted | 24.49 | 27.05 | 29.75 | 1.72 | 3.12 | 0.32
Table 1: Mean $l_{2}$ and Wasserstein distances between original and
adversarial images over 90 randomly selected images. Figure 2: Targeted
$l_{2}$ attacks (left) and Wasserstein attacks (right). Rows represent the
original while columns the target class.

For the numerical experiments, we use the notation of $l_{2}$ and Wasserstein
attacks based on the objective function in the model (1). For MNIST, we would
like to stress that none of the images was manually selected, and all images
were generated randomly.

Figure 2 shows targeted $l_{2}$ attacks (left) and targeted Wasserstein
attacks (right). Even though both attacks performed well, the Wasserstein
attack shows fewer grey artefacts around the digits. This is natural as the
Wasserstein distance considers the spatial distribution of pixels. The is not
the case for the $l_{2}$ distance.

Table 1 shows a numerical comparison between these two attacks and the CW
$l_{2}$ attack (Carlini and Wagner 2017) implemented in the Foolbox library
(Rauber et al. 2020). Besides the targeted attacks, we also implemented
untargeted attacks and attacks against the robust network (additional figures
are in the appendix). We show the $l_{2}$ and the Wasserstein distances. The
table shows that if we optimize with respect to the Wasserstein distance, the
Wasserstein distance between the original and adversarial images (column 8) is
the smallest. The same holds for the $l_{2}$ distance (column 4). Even though
the CW attack generated smaller values in the $l_{2}$ distance, this happened
because it is not restricted by the decoder. Moreover, as we will see from the
next figure, the CW attack generates lower-quality images which do not keep
the hidden structure of the original images. It is also not surprising that
the needed perturbations are smaller for untargeted attacks and the non-robust
network.

Figure 3: Comparison of CW and our attacks. Least significant bit (LSB) is
defined in (6), the difference is between the original and adversarial images.
The figure shows that out attacks keep the hidden structure of the data (LSB)
and that the perturbations are not located randomly but around edges.

Figure 3 shows a visual comparison for untargeted attacks. It shows the
original image (left) and the CW (middle left), $l_{2}$ (middle right) and
Wasserstein (right) attacks. Each three columns contain the adversarial image,
the least significant bit representation and the difference between the
original and adversarial images. Least significant bit is a standard
steganographic defence technique which captures changes in the hidden
structure of the data. For pixel intensities in $x\in[0,1]$, it is defined by

$\operatorname{lsb}(x)=\operatorname{mod}(\operatorname{round}(255x),2).$ (6)

Therefore, it transforms the image $x\in[0,1]$ into its $8$-bit representation
and takes the last (least significant) bit.

A significant difference between our attacks and the pixel-based CW attack is
that our attacks provide interpretability of perturbations. Figure 3 shows
that there is no clear pattern in the pixel perturbations of the CW attack,
neither in the least significant bit nor in the (almost) uniformly distributed
perturbations. On the other hand, our method keeps the least significant bit
structure. This implies that the least significant bit defence can easily
recognize the CW attack, while our attacks cannot be recognized. At the same
time, our methods concentrate the attack around the edges of the image.
Therefore, our attacks are compatible with the basic steganographic rule
stating that attacks should be concentrated mainly around key features such as
edges. We show the same figure for the robust classifier in the appendix.

### Numerical results for Imagenet

When fooling the ImageNet classifier, we use only the targeted version of our
algorithm to prevent trivial class changes in the case of the untargeted
attack, such as changing the dog’s breed. We manually select several pictures
of various animals as original images. As the target class we chose broccoli
due to its distinct features.

Figure 5 demonstrates the differences between original and adversarial images.
In all experiments, we perturbed the second intermediate layer of the BigBiGAN
decoder. The first row shows the original images, while the second row shows
the adversarial images misclassified all as broccoli. The third row highlights
the difference between original and adversarial images by taking the mean
square root of the absolute pixel difference across all channels. This
difference shows interpretable silhouettes, which allows us to assign
semantical meaning to many perturbations. Most of the differences are
concentrated in key components of the animals, such as changes of texture,
colouring and brightness of their semantically meaningful components. For
example, the dog fur has a different texture, and the grass has lighter
colour. Another example is the bird’s beak, which is longer in the third and
thicker in the fourth adversarial image. The difference plot also shows that
significant perturbations happen around animals edges or key features such as
eyes and nose.

Figure 4: Effect of the intermediate layer choice. The columns show the
perturbation in the first (left), second (middle) and third (right)
intermediate layer. Figure 5: Unperturbed image (top), perturbed image
(middle) and their difference (bottom). The perturbations happen around key
animal features. Figure 6: Development of perturbed images over iterations.
Broccoli features get gradually incorporated into the dog.

The crucial idea of our paper is to perturb the intermediate decoder layers.
However, we can place the perturbation in any intermediate layer, each with a
different effect. The earlier layers generate high-level features, while later
layers generate much finer features or even perturbations in pixels without
any semantic meaning. Figure 4 shows the impact of the choice of the
intermediate layer. The columns show the adversarial images when perturbing
the first (left column), second (middle column) and third (right column)
intermediate layer. The figure shows that the deeper we perturb the decoder,
the more the perturbations shift from high-level to finer perturbations. The
first layer results in a purely black bird and uniform background. The
silhouettes show that while the first layer perturbs the colour of the chest,
the beak or the background, the silhouette in the third layer disappears, and
the perturbations amount almost to pixel perturbations. We present additional
results of perturbing other layers in the appendix.

Figure 6 documents a similar effect of high-level perturbations. While the
previous figure showed the dependence of high-level features on the
intermediate layer choice, this figure shows this dependence on the
development of iterations in Algorithm 1. The left image shows the broccoli
image used for the initialization of our algorithm. The right image shows the
unperturbed dog. The middle five images visualize how the perturbations modify
the broccoli image as the number of iterations increases. At the beginning of
the algorithm run, the perturbations are interpretable: the dog is
reconstructed by high-level broccoli features. As the algorithm gradually
converges, the perturbed image increasingly resembles the original dog as the
broccoli features fuse with the dog image. Some of the original broccoli
features remain integrated into the final adversarial image. For example, the
next-to-last picture contains a slightly lighter green shade than the original
dog picture. Similarly, we can interpret the changes in the dog fur and muzzle
textures as they originate from the broccoli texture. We would like to point
out that all central images are classified as broccoli due to the construction
of our algorithm. We show the impact of the choice of the initial broccoli
image on the final adversarial image in the appendix.

##  5 Conclusion

This paper presented a novel method for generating adversarial images. Instead
of the standard pixel perturbation, we use an encoder-decoder model and
perturb high-level features from intermediate decoder layers. We showed that
our method generates high-quality adversarial images in both targeted and
untargeted settings on the MNIST and ImageNet datasets. We demonstrated that
our adversarial images indeed perturb high-level features. This makes them
more resilient to being recognized as adversarial by standard defence
techniques.

## References

  * Brendel, Rauber, and Bethge (2018) Brendel, W.; Rauber, J.; and Bethge, M. 2018.  Decision-Based Adversarial Attacks: Reliable Attacks Against Black-Box Machine Learning Models.  In _International Conference on Learning Representations (ICLR) 2018_.
  * Carlini and Wagner (2017) Carlini, N.; and Wagner, D. 2017.  Towards Evaluating the Robustness of Neural Networks.  In _2017 IEEE Symposium on Security and Privacy (SP)_ , 39–57.
  * Chen, Jordan, and Wainwright (2020) Chen, J.; Jordan, M. I.; and Wainwright, M. J. 2020.  HopSkipJumpAttack: A Query-Efficient Decision-Based Attack.  In _2020 IEEE Symposium on Security and Privacy (SP)_ , 1277–1294.
  * Cuturi (2013) Cuturi, M. 2013.  Sinkhorn Distances: Lightspeed Computation of Optimal Transport.  In _Advances in Neural Information Processing Systems 26_ , volume 26, 2292–2300.
  * Deng et al. (2009) Deng, J.; Dong, W.; Socher, R.; Li, L.-J.; Li, K.; and Fei-Fei, L. 2009.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE conference on computer vision and pattern recognition_ , 248–255.
  * Donahue, Krähenbühl, and Darrell (2016) Donahue, J.; Krähenbühl, P.; and Darrell, T. 2016.  Adversarial feature learning.  _arXiv preprint arXiv:1605.09782_ .
  * Donahue and Simonyan (2019) Donahue, J.; and Simonyan, K. 2019.  Large Scale Adversarial Representation Learning.  In _Advances in Neural Information Processing Systems_ , volume 32, 10541–10551.
  * Dumoulin et al. (2016) Dumoulin, V.; Belghazi, I.; Poole, B.; Mastropietro, O.; Lamb, A.; Arjovsky, M.; and Courville, A. 2016.  Adversarially learned inference.  _arXiv preprint arXiv:1606.00704_ .
  * Feydy et al. (2019) Feydy, J.; Séjourné, T.; Vialard, F.-X.; Amari, S.-i.; Trouve, A.; and Peyré, G. 2019.  Interpolating between Optimal Transport and MMD using Sinkhorn Divergences.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2681–2690.
  * Goodfellow, Shlens, and Szegedy (2015) Goodfellow, I. J.; Shlens, J.; and Szegedy, C. 2015.  Explaining and Harnessing Adversarial Examples.  In _International Conference on Learning Representations (ICLR) 2015_.
  * Johnson and Jajodia (1998) Johnson, N. F.; and Jajodia, S. 1998.  Exploring steganography: Seeing the unseen.  _Computer_ 31(2): 26–34.
  * Kurakin, Goodfellow, and Bengio (2016) Kurakin, A.; Goodfellow, I.; and Bengio, S. 2016.  Adversarial machine learning at scale.  _arXiv preprint arXiv:1611.01236_ .
  * LeCun, Cortes, and Burges (2010) LeCun, Y.; Cortes, C.; and Burges, C. 2010.  MNIST handwritten digit database.  _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ 2\.
  * Madry et al. (2018) Madry, A.; Makelov, A.; Schmidt, L.; Tsipras, D.; and Vladu, A. 2018.  Towards Deep Learning Models Resistant to Adversarial Attacks.  In _International Conference on Learning Representations (ICLR) 2018_.
  * Meng and Chen (2017) Meng, D.; and Chen, H. 2017.  MagNet: A Two-Pronged Defense against Adversarial Examples.  In _Proceedings of the 2017 ACM SIGSAC Conference on Computer and Communications Security_ , 135–147.
  * Nocedal and Wright (2006) Nocedal, J.; and Wright, S. 2006.  _Numerical optimization_.  Springer Science & Business Media.
  * Papernot et al. (2016) Papernot, N.; McDaniel, P.; Jha, S.; Fredrikson, M.; Celik, Z. B.; and Swami, A. 2016.  The Limitations of Deep Learning in Adversarial Settings.  In _2016 IEEE European Symposium on Security and Privacy (EuroS &P)_, 372–387.
  * Rauber et al. (2020) Rauber, J.; Zimmermann, R.; Bethge, M.; and Brendel, W. 2020.  Foolbox Native: Fast adversarial attacks to benchmark the robustness of machine learning models in PyTorch, TensorFlow, and JAX.  _Journal of Open Source Software_ 5(53): 2607.  doi:10.21105/joss.02607.  URL https://doi.org/10.21105/joss.02607.
  * Samangouei, Kabkab, and Chellappa (2018) Samangouei, P.; Kabkab, M.; and Chellappa, R. 2018.  Defense-GAN: Protecting Classifiers Against Adversarial Attacks Using Generative Models.  In _International Conference on Learning Representations (ICLR) 2018_.
  * Simonyan and Zisserman (2015) Simonyan, K.; and Zisserman, A. 2015.  Very Deep Convolutional Networks for Large-Scale Image Recognition.  In _International Conference on Learning Representations (ICLR) 2015_.
  * Song et al. (2018) Song, Y.; Shu, R.; Kushman, N.; and Ermon, S. 2018.  Constructing Unrestricted Adversarial Examples with Generative Models.  In _Advances in Neural Information Processing Systems_ , volume 31, 8312–8323.
  * Szegedy et al. (2014) Szegedy, C.; Zaremba, W.; Sutskever, I.; Bruna, J.; Erhan, D.; Goodfellow, I.; and Fergus, R. 2014.  Intriguing properties of neural networks.  In _International Conference on Learning Representations (ICLR) 2014_.
  * Tan and Le (2019) Tan, M.; and Le, Q. V. 2019.  EfficientNet: Rethinking Model Scaling for Convolutional Neural Networks.  In _International Conference on Machine Learning_ , 6105–6114.
  * Wong, Schmidt, and Kolter (2019) Wong, E.; Schmidt, F. R.; and Kolter, J. Z. 2019.  Wasserstein Adversarial Examples via Projected Sinkhorn Iterations.  In _International Conference on Machine Learning_ , 6808–6817.
  * Xiao et al. (2018) Xiao, C.; Li, B.; yan Zhu, J.; He, W.; Liu, M.; and Song, D. 2018.  Generating Adversarial Examples with Adversarial Networks.  In _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence_ , 3905–3911.

The Appendix presents the proof of Theorem 1 and shows additional results.

##  Appendix A Proof of Theorem 1

Recall that a function $f$ is Lipschitz continuous with constant $L$ if

$\|f(x)-f(y)\|\leq L\|x-y\|$

for all $x$ and $y$.

######  Theorem 2.

Assume that $g\circ D_{2}$ is a continuous function. Then Algorithm 2
generates a strictly feasible point $p$ in a finite number of iterations. If
$g\circ D_{2}$ is moreover Lipschitz continuous with modulus $L$, then
Algorithm 2 converges in at most $\log_{2}\frac{L\|p\|}{\delta}$ iterations.

###### Proof.

Denote the iterations from Algorithm 2 by $b^{k}$ and $a^{k}$ with the
initialization $b^{0}=1$ and $a^{0}=0$. Since the interval halves at each
iteration, we obtain

$b^{k}-a^{k}=\frac{1}{2}(b^{k-1}-a^{k-1})=\frac{1}{2^{k}}(b^{0}-a^{0})=\frac{1}{2^{k}}.$
(7)

Define

$\tilde{g}(c)=\hat{g}((1-c)p+cp_{\rm next}).$

Since we have

$\displaystyle\tilde{g}(b^{0})$ $\displaystyle=\tilde{g}(1)=\hat{g}(p_{\rm
next})\leq 0,$ $\displaystyle\tilde{g}(a^{0})$
$\displaystyle=\tilde{g}(0)=\hat{g}(p)<0.$

due to the construction of the algorithm, we have $\tilde{g}(b^{k})\leq 0$ and
$\tilde{g}(a^{k})>0$ for all $k$. Since $\tilde{g}$ is a continuous function
due to continuity of $D_{2}\circ g$, Algorithm 2 converges in a finite number
of iterations.

Assume that $g\circ D_{2}$ is a Lipschitz continuous function with modulus
$L$, then $\hat{g}$ is also Lipschitz continuous with modulus $L$. Then

$\displaystyle|\tilde{g}$ $\displaystyle(c_{1})-\tilde{g}(c_{2})|$
$\displaystyle=|\hat{g}((1-c_{1})p+c_{1}p_{\rm
next})-\hat{g}((1-c_{2})p+c_{2}p_{\rm next})|$ $\displaystyle\leq
L\|(1-c_{1})p+c_{1}p_{\rm next}-(1-c_{2})p-c_{2}p_{\rm next}\|$
$\displaystyle=L\|(p_{\rm next}-p)(c_{1}-c_{2})\|$ $\displaystyle\leq
L\|p_{\rm next}-p\||c_{1}-c_{2}|$

Therefore, $\tilde{g}$ is a Lipschitz continuous function with constant
$L\|p_{\rm next}-p\|$. Together with (7) this implies

$\displaystyle\tilde{g}(b^{k})-\tilde{g}(a^{k})$ $\displaystyle\leq L\|p_{\rm
next}-p\|(b^{k}-a^{k})$ (8) $\displaystyle=L\|p_{\rm
next}-p\|\frac{1}{2^{k}}.$

If the algorithm did not finish within $k$ iterations, we have

$\tilde{g}(b^{k})-\tilde{g}(a^{k})\geq\delta.$ (9)

The comparison of (8) and (9) implies that the algorithm cannot run for more
than $\log_{2}\frac{L\|p_{\rm next}-p\|}{\delta}$ iterations, which finishes
the proof. ∎

##  Appendix B Additional results

This section extends the results from the main manuscript body. We will always
present a figure and then compare it with the corresponding figure from the
manuscript body. The former figures start with a letter while the latter
figures with a digit.

Figure 7 shows the untargeted attacks for the non-robust network. It
corresponds to Figure 2 from the manuscript body. The images are again nice,
with the Wasserstein attack performing better than the $l_{2}$ attack. The
small digit in each subfigure shows to which class the digit was
misclassified. As we have already mentioned, our method always works with
feasible points and, therefore, all digits were successfully misclassified. In
other words, these images were generated randomly without the need to select
them manually.

Figure 8 shows the attacks on the robust network of (Madry et al. 2018). It
corresponds to Figure 3 from the main manuscript body. The quality of the CW
attack increased. It now preserves the least significant bit structure, and
the perturbations shifted towards the digit edges. The Wasserstein attack
keeps the superb performance with visually the same results as in Figure 3. We
conclude that our attacks are efficient against adversarially trained
networks.

Figure 7: Untargeted $l_{2}$ attacks (left) and Wasserstein attacks (right).
Rows represent the class. The small number in each subfigure shows to which
class the digit was classified.

Figure 9 shows the effect of the choices of the intermediate layer to perturb
and of the initial point for Algorithm 1. It corresponds to Figure 4 from the
main manuscript body. The first column shows the point which was used to
initialize Algorithm 1. The next three columns present the adversarial images
with the perturbation inserted in different intermediate layers. The last
column shows the original image. The effect of the initialization is
negligible when perturbing the third intermediate layer (column 4). However,
it has a huge impact when perturbing earlier intermediate layers. The
intermediate latent representation of the second broccoli (rows 2 and 4)
carries a preference for the white colour, which is visible in the adversarial
image for the first intermediate layer (column 2).

Figure 10 also shows the effect of the initial point for Algorithm 1. It
corresponds to Figure 6 from the main manuscript body. We see that the
features of the initial broccoli (column 1) gradually incorporate into the dog
image (columns 2-5). The adversarial image (column 6) still have some
connection to the initial broccoli, for example, in the background colour. The
adversarial image is close to the original image (column 7). This figure shows
once again that our algorithm works with high-level features and not pixel
modifications.

Figure 8: Comparison of CW and our attacks on the robust network. Least
significant bit (LSB) is defined in (6), the difference is between the
original and adversarial images. Figure 9: Effect of the intermediate layer
choice. The columns show initial point for Algorithm 1 (column 1), the
adversarial image when perturbations were performed in the first (column 2),
second (column 3) and third (column 4) intermediate layer and the original
image (column 5).

Figure 10: Development of perturbed images over iterations. Broccoli features
get gradually incorporated into the dog.