]Faculty of Engineering and Physical Sciences, University of Southampton,
United Kingdom

# Sensitivity of long periodic orbits of chaotic systems

D. Lasagna  [

###### Abstract

The properties of long, numerically-determined periodic orbits of two low-
dimensional chaotic systems, the Lorenz equations and the Kuramoto-Sivashinsky
system in a minimal-domain configuration, are examined. The primary question
is to establish whether the sensitivity of period averaged quantities with
respect to parameter perturbations computed over long orbits can be used as a
sufficiently good proxy for the response of the chaotic state to finite-
amplitude parameter perturbations. To address this question, an inventory of
thousands of orbits at least two-order of magnitude longer than the shortest
admissible cycles is constructed. The expectation of period averages, Floquet
exponents and sensitivities over such set is then obtained. It is shown that
all these quantities converge to a limiting value as the orbit period is
increased. However, while period averages and Floquet exponents appear to
converge to analogous quantities computed from chaotic trajectories, the
limiting value of the sensitivity is not necessarily consistent with the
response of the chaotic state, similar to observations made with other
shadowing algorithms.

Unstable Periodic Orbits, Adjoint Sensitivity, Floquet analysis

###### pacs:

Valid PACS appear here

††preprint: AIP/123-QED

##  I Introduction

Evidence has been offered in recent years Viswanath (2007); Chandler and
Kerswell (2013); Willis  _et al._ (2016); Budanur  _et al._ (2017) that
temporally-recurrent invariant solutions of the Navier-Stokes equations –
unstable periodic orbits – may provide a constitutive skeleton organizing
spatiotemporal dynamics of turbulent shear flows in canonical geometries.
Motivated by these advances, we have recently suggested Lasagna (2018) how
unstable periodic orbits may be used to design control strategies for shear
flows, rather than serving as a tool to rationalise turbulence dynamics. In
particular, we have specialized adjoint methods for time-periodic systems
Hwang and Choi (2008); Giannetti  _et al._ (2010); Meliga  _et al._ (2016);
Giannetti  _et al._ (2019) to unstable periodic orbits. We have shown that
enforcing periodicity conditions on the adjoint problem, justified by the
peculiar topology of these trajectories, prevents the growth of exponential
instabilities that would otherwise feature prominently in the solution of the
adjoint equations Lea  _et al._ (2000); Luchini and Quadrio (2002).

Operationally, the approach provides the sensitivity of period averaged
quantities with respect to small perturbations of variables parametrising the
equations of motion. Geometrically, small parameter perturbations can be
pictured as producing smooth, global state-space deformations of the unstable
periodic orbits supporting and shaping the attractor, as opposed to causing
exponential divergence. The approach is a special case of shadowing theory
ideas Bowen (1975); Ruelle (1999); Nusse and Yorke (1988), recently introduced
in the context of sensitivity analysis of chaotic systems Wang (2013); Wang
_et al._ (2014); Ni and Wang (2017); Shawki and Papadakis (2018).

In principle, complete knowledge of the short, fundamental cycles should
suffice to compute ergodic averages using cycle averaging formulae Eckhardt
and Ott (1994); Cvitanović (1995); Artuso  _et al._ (1999). Of relevance for
our original motivation is that a formalism that relies on the sensitivity of
such cycles to compute the sensitivity of ergodic averages was proposed in
Refs. Kazantsev (1998, 2001). Obtaining all short cycles up to a given
topological length may be practical for low-dimensional systems (see e.g. Ref.
Viswanath (2003)). However, this step has proven more challenging for
turbulent shear flows Chandler and Kerswell (2013); Lucas and Kerswell (2015),
given well documented difficulties in locating invariant solutions van Veen
_et al._ (2019). This issue is particularly relevant, since the quality of
cycle averaging predictions using incomplete hierarchies is as good as the
most important orbit that one fails to locate Budanur  _et al._ (2015).

In light of such issues, we explore in this paper a heuristic approach whereby
available computational resources are spent to locate one or few periodic
orbits, with sufficiently long period for them to span a good fraction of the
attractor. The open question is whether sensitivities of period averaged
quantities computed over a long orbit can be useful approximations for the
response of ergodic averages to parameter perturbations. Like period averages,
sensitivities computed over different periodic orbits vary from cycle to
cycle. However, it is known that the period average of long periodic orbits
converges to a defined value when the period increase Zaks and Goldobin
(2010). Evidence showing a similar convergence for sensitivities of period
averages is currently not available in the literature and would provide
initial support for the above heuristic. This approach is, admittedly, guided
more by empiricism rather than by a solid theoretical basis. Hence, the aim of
this paper is to make a first step in exploring its viability. We resort, by
necessity, to low-dimensional systems, where obtaining a sizeable inventory of
long periodic orbits is feasible. We consider the Lorenz equations at standard
parameters Lorenz (1963) and a small-domain Kuramoto-Sivashinsky system in the
antisymmetric subspace Christiansen  _et al._ (1997); Lan and Cvitanović
(2008).

One remark is in order. It is true that some long periodic orbits might not
provide good approximations, for instance orbits close to bifurcation or
orbits visiting certain areas of the attractor where the response is
particularly large. As an illustrative example, in Ref. Lasagna (2018), long
periodic orbits of the Lorenz equations passing close to the origin were found
that have extremal value of the sensitivity. Hence, to develop a quantitative
understanding of the potential impact of such extremal orbits, we locate
thousands of long periodic orbits by converging near recurrences explored by
the chaotic flow and report sensitivity predictions of these periodic orbits
in a statistical manner. The analysis of this ensemble is clearly bound to
suffer of some form of survivor bias. For instance, Newton-Raphson techniques
often fail to converge when the orbit shadowed by the near recurrence event
possesses multipliers close to the unit circle Sánchez  _et al._ (2004).
Nevertheless, this bias does not prevent us to complete the original task,
which is to analyse the properties of typical periodic orbits that may be
found numerically rather than the properties of all admissible cycles.

The rest of this paper is organized as follows. In section II we lay down the
general notation. This is followed by section III, where we recall fundamental
elements of Floquet theory for the linear stability of periodic orbits Farkas
(1994); Moore (2006). In section IV we recall the sensitivity technique for
unstable periodic orbits originally introduced in Ref. [Lasagna, 2018]. With
this introductory material, in section V we establish the connection between
Floquet stability and sensitivity of periodic orbits. This analysis
illustrates more clearly the inevitable effects of bifurcations that periodic
orbits undergo as parameters are varied and supports the interpretation of
numerical results, reported in section VI.

##  II Preliminaries

We consider dissipative dynamical systems of the form

$\frac{\mathrm{d}\mathbf{u}(t)}{\mathrm{d}t}=\mathbf{f}(\mathbf{u}(t),\alpha),$
(1)

governing the evolution of the state vector $\mathbf{u}(t)\in\mathbb{R}^{n}$,
with $t$ being time. We restrict the attention to problems where the nonlinear
vector field $\mathbf{f}$ depends only on one parameter $\alpha\in\mathbb{R}$,
in regimes for which chaotic solutions are observed for typical initial
conditions. For systems that depend on multiple parameters, the sensitivity of
time averaged quantities can be analysed over each parameter independently, or
with an adjoint approach. We denote with
$\mathbf{f}_{\mathbf{u}}(t)\in\mathbb{R}^{n\times n}$ the stability matrix

$\mathbf{f}_{\mathbf{u}}(t)=\mathbf{f}_{\mathbf{u}}(\mathbf{u}(t),\alpha)=\frac{\partial\mathbf{f}(\mathbf{u}(t),\alpha)}{\partial\mathbf{u}}.$
(2)

and define the Jacobian matrices $\mathbf{M}(t,\tau)\in\mathbb{R}^{n\times
n}$, $t\geq\tau$, satisfying the initial value problem

$\frac{\mathrm{d}\mathbf{M}(t,\tau)}{\mathrm{d}t}=\mathbf{f}_{\mathbf{u}}(t)\mathbf{M}(t,\tau),\quad\mathbf{M}(\tau,\tau)=\mathbf{I},$
(3)

with $\mathbf{I}$ the identity matrix. The dot product of two vectors is
denoted with $\mathbf{a}^{\top}\cdot\mathbf{b}$, with the script ⊤ indicating
transposition of vectors and matrices.

We focus on periodic solutions of (1), satisfying
$\mathbf{u}(t+T)=\mathbf{u}(t)$ for an unknown period $T$ that is not set a
priori, but depends implicitly on the parameter $\alpha$. We shall thus
consider the space of smooth periodic functions

$\mathcal{P}_{T}=\\{f(t):\mathbb{R}\mapsto\mathbb{R},f(t)=f(t+T)\\},$ (4)

parametrized by $T$, and extend this space to vector-valued functions, denoted
with $\mathcal{P}_{T}^{n}$. We will make use of the norm
$\|\cdot\|_{\mathcal{P}_{T}^{n}}$ induced by the inner product

$\left\llbracket\mathbf{v}(t),\mathbf{w}(t)\right\rrbracket=\frac{1}{T}\int_{0}^{T}\mathbf{v}(t)^{\top}\cdot\mathbf{w}(t)\,\mathrm{d}t,$
(5)

for any two vector-valued functions $\mathbf{v}(t)$ and $\mathbf{w}(t)$ in
$\mathcal{P}_{T}^{n}$.

##  III Elements of Floquet theory

In this section, we briefly recall some elements of Floquet theory
Guckenheimer and Holmes ; Moore (2006) and define the direct and adjoint
Floquet eigenfunctions. These functions are used in section III as invariant
subspaces in which the solution of the sensitivity problem can be conveniently
expanded, shedding light on the relation between stability and sensitivity of
period averages.

For an arbitrary point on a periodic orbit, the eigenvalue decomposition of
the monodromy matrix

$\mathbf{M}(T,0)\mathbf{e}_{k}(0)=\mu_{k}\mathbf{e}_{k}(0),\quad
k=1,\ldots,n,$ (6)

produces the Floquet multipliers $\mu_{k}$ and the associated right
eigenvectors $\mathbf{e}_{k}(0)$. To simplify the notation and the analysis of
section V, we consider here the special case in which multipliers are real and
positive. In more general cases, new function spaces in addition to
$\mathcal{P}_{T}^{n}$ are required (see Ref. Moore (2006)). Assuming the
multipliers are distinct, the eigenvectors form collectively a basis of
$\mathbb{R}^{n}$. From the multipliers, the Floquet exponents
$\lambda_{k}=\log(\mu_{k})/T$ can be calculated, defining the period averaged
growth/decay of tangent perturbations initially aligned to the invariant
subspaces $\mathbf{e}_{k}(0)$.

We introduce the Floquet eigenfunctions
$\mathbf{w}_{k}(t)\in\mathcal{P}_{T}^{n}$, generated by the eigenvectors
$\mathbf{e}_{k}(0)$ as

$\mathbf{w}_{k}(t)=\exp\left(-\lambda_{k}t\right)\mathbf{M}(t,0)\mathbf{e}_{k}(0),$
(7)

and satisfying the differential eigenvalue problem

$\mathcal{L}\mathbf{w}_{k}(t)\equiv\frac{\mathrm{d}\mathbf{w}_{k}(t)}{\mathrm{d}t}-\mathbf{f}_{\mathbf{u}}(t)\mathbf{w}_{k}(t)=-\lambda_{k}\mathbf{w}_{k}(t),$
(8)

with unit $\|\cdot\|_{\mathcal{P}_{T}^{n}}$ norm. Sorting the Floquet
exponents in descending order, we denote by $\chi-1$ the number of positive
exponents. It is well known that
$\mathbf{w}_{\chi}(t)=\mathbf{f}(t)/{\|\mathbf{f}(t)\|_{\mathcal{P}_{T}^{n}}}$
is a marginal direction. In other words, the linear differential operator
$\mathcal{L}$ is singular, with nullspace

$\mathrm{Null}\left\\{\mathcal{L}\right\\}=\mathrm{span}\\{{\mathbf{w}_{\chi}}(t)\\}.$
(9)

The other useful element of Floquet theory for our purposes, perhaps
considered less extensively in the literature Marcotte and Grigoriev (2016),
are the adjoint Floquet eigenfunctions $\mathbf{w}^{+}_{k}(t)$. These are
elements of $\mathcal{P}_{T}^{n}$ defined as

$\mathbf{w}^{+}_{k}(t)=\exp(-\lambda_{k}(T-t))\mathbf{M}(t,T)^{\top}\mathbf{e}_{k}^{+}(T),$
(10)

where the vectors $\mathbf{e}^{+}_{k}(T)$ are the left eigenvectors of the
monodromy matrix, satisfying

$\mathbf{e}^{+}_{k}(T)^{\top}\mathbf{M}(0,T)=\mu_{k}\mathbf{e}^{+}_{k}(T)^{\top}$
(11)

for the same multipliers and exponents of the direct problem. The adjoint
eigenfunctions satisfy the adjoint differential eigenproblem

$\mathcal{L}^{+}\mathbf{w}_{k}^{+}(t)\equiv-\frac{\mathrm{d}\mathbf{w}_{k}^{+}(t)}{\mathrm{d}t}-\mathbf{f}_{\mathbf{u}}^{\top}(t)\mathbf{w}_{k}^{+}(t)=-\lambda_{k}\mathbf{w}^{+}_{k}(t),$
(12)

where the operator $\mathcal{L}^{+}$ is the adjoint of $\mathcal{L}$ according
to the inner product (5). These two operators share the same spectrum and the
adjoint operator $\mathcal{L}^{+}$ is thus singular with nullspace

$\mathrm{Null}\left\\{\mathcal{L}^{+}\right\\}=\mathrm{span}\\{{\mathbf{w}^{+}_{\chi}}(t)\\}.$
(13)

A final remark is that, at any time $t$, the direct or adjoint eigenfunctions
do not form individually an orthogonal set of vectors but instead satisfy the
bi-orthogonality relation

$\mathbf{w}_{k}^{+}(t)^{\top}\cdot\mathbf{w}_{j}(t)=\delta_{kj}C_{k}\quad\forall
t,$ (14)

with $\delta_{kj}$ the Kronecker symbol and with $C_{k}\neq 0\in\mathbb{R}$.

##  IV Sensitivity analysis

We now briefly recall the tangent approach to compute the sensitivity of
period averaged quantities Lasagna (2018). Consider an observable of interest,
denoted by a function

$J(t)=J(\mathbf{u}(t)):\mathbb{R}^{n}\mapsto\mathbb{R},$ (15)

that, for the sake of simplicity, is assumed here not to depend explicitly on
the parameter $\alpha$. Let also denote the partial derivative of the
observable with respect to its argument as
$J_{\mathbf{u}}(t)=J_{\mathbf{u}}(\mathbf{u}(t))$, mapping $\mathbb{R}^{n}$ to
$\mathbb{R}^{n}$. Consider one periodic orbit
$\mathbf{u}(t)\in\mathcal{P}_{T}^{n}$ and define the function
$\mathcal{J}(\alpha):\mathbb{R}\mapsto\mathbb{R}$

$\displaystyle\mathcal{J}(\alpha)=\frac{1}{T}\int_{0}^{T}J(\mathbf{u}(t))\mathrm{d}t$
(16)

as the period average of the observable over the periodic orbit. This is a
function of the parameter $\alpha$ since both $\mathbf{u}(t)$ and $T$ depend,
implicitly, on it. The goal is to compute the sensitivity of the period
average with respect to $\alpha$, the gradient
${\mathrm{d}\mathcal{J}}/{\mathrm{d}\alpha}$ (also shortened to
$\mathcal{J}_{\alpha}$ in the following sections).

Linearisation of (16) as reported in Ref. [Lasagna  _et al._ , 2019] shows
that the gradient of the period average is given by the inner product

${\mathrm{d}\mathcal{J}}/{\mathrm{d}\alpha}=\left\llbracket
J_{\mathbf{u}}(t),\mathbf{v}(t)\right\rrbracket,$ (17)

where the perturbation $\mathbf{v}(t)\in\mathcal{P}_{T}^{n}$ satisfies the
tangent equation

$\mathcal{L}\mathbf{v}(t)=\mathbf{f}_{\alpha}(t)-\tau\mathbf{f}(t).$ (18)

The perturbation $\mathbf{v}(t)$ is the first order state-space deformation of
the periodic orbit $\mathbf{u}(t)$ when $\alpha$ is varied. The scalar
$\tau=(\mathrm{d}T/\mathrm{d}\alpha)/T$ is the (unknown) relative period
gradient, producing an algebraically growing mode along $\mathbf{f}(t)$,
allowing $\mathbf{v}(t)$ to be time periodic. The introduction of this term is
akin to classical approaches in perturbation/continuation analysis of periodic
problems Farkas (1994); Viswanath (2001), where time is rescaled by the
period. The forcing term $\mathbf{f}_{\alpha}(t)\in\mathcal{P}^{n}_{T}$ is the
derivative of the right hand side of (1) with respect to the parameter

$\mathbf{f}_{\alpha}(t)=\mathbf{f}_{\alpha}(\mathbf{u}(t),\alpha)=\frac{\partial\mathbf{f}(\mathbf{u}(t),\alpha)}{\partial\alpha}.$
(19)

For an hyperbolic orbit, the differential operator $\mathcal{L}$ is singular
with nullspace given by (9) and equation (18) has a one parameter family of
solutions. Physically speaking, this is a reflection of the translational
invariance along a periodic orbit: if $\mathbf{v}(t)$ is a solution, then
$\mathbf{v}(t)+\sigma\mathbf{f}(t)$ is a solution too, for any
$\sigma\in\mathbb{R}$, with same gradient
$\mathrm{d}\mathcal{J}/\mathrm{d}{\alpha}$, since

$\llbracket J_{\mathbf{u}}(t),\mathbf{f}(t)\rrbracket=0$ (20)

for all possible cost functions when $\mathbf{u}(t)\in\mathcal{P}^{n}_{T}$.
Hence, for (18) to have a solution, the scalar $\tau$ must have the unique
value that shifts the right hand side
$\mathbf{f}_{\alpha}(t)-\tau\mathbf{f}(t)$ in the range of the operator
$\mathcal{L}$, or, by Fredholm’s Alternative (cf. [Hale, 1969], Lemma 1.1, pg.
146) makes it orthogonal to the nullspace of its adjoint $\mathcal{L}^{+}$,
i.e. by satisfying

$\left\llbracket\mathbf{w}_{\chi}^{+}(t),\mathbf{f}_{\alpha}(t)-\tau\mathbf{f}(t)\right\rrbracket=0.$
(21)

Hence, the scalar $\tau$ could be in principle determined as

$\tau=\frac{\left\llbracket\mathbf{w}_{\chi}^{+}(t),\mathbf{f}_{\alpha}(t)\right\rrbracket}{C_{\chi}\|\mathbf{f}(t)\|_{\mathcal{P}^{n}_{T}}}$
(22)

if $\mathbf{w}_{\chi}^{+}(t)$ was known. Numerically, it is more convenient to
drop the singularity by adding the constraint

$\left\llbracket\mathbf{v}(t),\mathbf{f}(t)\right\rrbracket=0,$ (23)

which fixes the component of $\mathbf{v}(t)$ along the nullspace and leads to
the solution of (18) with minimum norm Moore (2006). In matrix form, the
tangent problem reads

$\displaystyle\left[\begin{array}[]{c|c}\displaystyle\mathcal{L}&\mathbf{f}(t)\\\
\hline\cr\\!\rule{0.0pt}{9.95863pt}\left\llbracket\cdot,\mathbf{f}(t)\right\rrbracket&0\\\
\end{array}\right]\cdot\left[\begin{array}[]{c}\mathbf{v}(t)\\\
\hline\cr\\!\rule{0.0pt}{9.95863pt}\tau\end{array}\right]=\left[\begin{array}[]{c}\mathbf{f}_{\alpha}(t)\\\
\hline\cr\\!\rule{0.0pt}{9.95863pt}0\end{array}\right]$ (24)

and its solution provides the perturbation $\mathbf{v}(t)$ and the period
gradient $\tau$. Note that the left hand side of (24) has the same structure
of the Newton-Raphson linear problems arising in the search of periodic
solutions Lan and Cvitanović (2004); Lasagna (2018), and similar
discretization techniques can be employed.

Constraining $\mathbf{v}(t)$ to remain in $\mathcal{P}^{n}_{T}$ by using an
appropriate numerical method Lasagna (2018) is the key to avoid exponential
instabilities intrinsic to the tangent dynamics around an unstable periodic
orbit Lea  _et al._ (2000). The solution $\mathbf{v}(t)$ will thus not grow
exponentially along the most unstable subspace $\mathbf{w}_{1}(t)$, but will
remain bounded, with magnitude and structure that depend on the complete
stability spectrum, as we shall see in section V. With a bounded
$\mathbf{v}(t)$, the gradient (17) is effectively the slope of the function
$\mathcal{J}(\alpha)$ obtained from continuation.

##  V Stability and sensitivity

In this section, we explain the relation between the linear stability of
periodic orbits and the sensitivity of period averages. We do not make use of
the following results for our numerical calculations in section VI, but aim to
develop tools to facilitate their interpretation.

Fundamentally, the approach consists in projecting the tangent problem (18)
onto the invariant subspaces formed by the Floquet eigenfunctions. With the
same assumptions on the multipliers as in section III, the solution
$\mathbf{v}(t)$ is expanded in the Floquet eigenfunctions,

$\mathbf{v}(t)=\sum_{k=1}^{n}\mathbf{w}_{k}(t)a_{k}(t),$ (25)

with unknown expansion coefficients $a_{k}(t)\in\mathcal{P}_{T}$. The forcing
term in (18) is also similarly expanded

$\mathbf{f}_{\alpha}(t)=\sum_{k=1}^{n}\mathbf{w}_{k}(t)b_{k}(t),$ (26)

where the functions $b_{k}(t)\in\mathcal{P}_{T}$ can be determined by dotting
(26) with the $k$-th adjoint eigenfunction

$b_{k}(t)=\frac{\mathbf{w}_{k}^{+}(t)^{\top}\cdot\mathbf{f}_{\alpha}(t)}{C_{k}},\quad
k=1,\ldots,n,$ (27)

where the bi-orthogonality relation (14) is used.

Substituting the expansion (25) into the tangent equation (18) and using (8),
produces

$\sum_{k=1}^{n}\mathbf{w}_{k}(t)\left[\frac{\mathrm{d}a_{k}(t)}{\mathrm{d}t}-\lambda_{k}a_{k}(t)-b_{k}(t)+\tau\|\mathbf{f}(t)\|_{\mathcal{P}^{n}_{T}}\delta_{k,\chi}\right]=\mathbf{0}.$
(28)

Since the Floquet eigenfunctions form a basis of $\mathbb{R}^{n}$ for all $t$
by assumption, the term in the square brackets in (28) must be zero. We thus
obtain a set of decoupled linear ODEs with constant coefficients

$\frac{\mathrm{d}a_{k}(t)}{\mathrm{d}t}=\lambda_{k}a_{k}(t)+b_{k}(t)-\tau\|\mathbf{f(t)}\|_{\mathcal{P}^{n}_{T}}\delta_{k,\chi},$
(29)

$k=1\ldots,n$, the tangent sensitivity problem expressed in the basis of the
Floquet eigenfunctions. Along the neutral subspace $\mathbf{f}(t)$ the
equation reads

$\frac{\mathrm{d}a_{\chi}(t)}{\mathrm{d}t}=\frac{\mathbf{w}_{\chi}^{+}(t)^{\top}\cdot\mathbf{f}_{\alpha}(t)}{C_{\chi}}-\tau\|\mathbf{f}(t)\|_{\mathcal{P}^{n}_{T}}.$
(30)

In order for $a_{\chi}(t)$ to remain in $\mathcal{P}_{T}$, the right hand side
must have zero integral over the period by Fredholm’s Alternative, since the
equation is self adjoint and $a_{\chi}^{+}(t)=1$ is a nontrivial solution of
the homogeneous adjoint equation. This constraints fixes the relative period
gradient $\tau$ to a value that is the same as equation (22).

A particularly insightful expression can be derived for the expanding and
contracting directions. The solution of the scalar ODEs
${\mathrm{d}a_{k}(t)}/{\mathrm{d}t}=\lambda_{k}a_{k}(t)+b_{k}(t)$ can be
expressed with a Green’s function approach (cf. [Hale, 1969], pg. 148) as

$a_{k}(t)=\int_{0}^{T}\frac{\mu_{k}}{1-\mu_{k}}e^{-\lambda_{k}s}b_{k}(t+s)\,\mathrm{d}s,$
(31)

and the upper bounds

$\displaystyle\sup_{t}|a_{k}(t)|\leq\;$
$\displaystyle\sup_{t}|b_{k}(t)|/{|\lambda_{k}|}=B_{k}/{|\lambda_{k}|},$ (32)

the key result of this section, can be derived.

This bound suggests several remarks. First, without further details on the
coefficients $b_{k}(t)$, generic parameter perturbations induce relatively
small state-space changes along the highly contracting or expanding
directions, while most of the “yield” occurs along the Floquet invariant
subspace associated to Floquet exponents with small magnitude. This is in
stark contrast with classical sensitivity analysis methods for chaotic
trajectories Lea  _et al._ (2000); Eyink  _et al._ (2004), where only the most
unstable covariant Lyapunov vector Ginelli  _et al._ (2007); Wolfe and
Samelson (2007) dominates asymptotically the solution of the tangent equations
Lea  _et al._ (2000). When varying the system parameter $\alpha$ towards a
bifurcation, one (or a pair of) Floquet exponent crosses the imaginary axis
and the tangent solution displays a large amplitude along the corresponding
direction, resulting in large gradients of time averaged quantities.

Second, the bound (32) shows that the amplitude of the tangent solution along
a particular Floquet eigenfunction $\mathbf{w}_{k}(t)$ depends directly on the
strength of the projection of the forcing $\mathbf{f}_{\alpha}(t)$ on the
associated adjoint Floquet eigenfunctions, the coefficients $b_{k}$. Hence,
for spatially extended systems, knowledge of the spatiotemporal dynamics of
the adjoint Floquet eigenfunctions and not just the direct ones Xu and Paul
(2018); Nikitin (2018); Inubushi  _et al._ (2015), might provide an
understanding of how physically relevant features of the solutions are
influenced by problem parameters Marcotte and Grigoriev (2016). One can then
interpret the adjoint Floquet eigenfunctions as special directions where the
forcing $\mathbf{f}_{\alpha}(t)$ can be particularly effective in modifying
dynamical behaviour Luchini and Bottaro (2014), which is useful, for instance,
for control design Giannetti  _et al._ (2019).

Third, and most importantly, the boundedness of the forcing term
$\mathbf{f}_{\alpha}(t)$ and of the Floquet eigenfunctions implies that the
coefficients $b_{k}(t)$, and thus the expansion coefficients $a_{k}(t)$ and
the tangent solution $\mathbf{v}(t)$, have, on average, similar magnitude for
long periodic orbits if the exponents of long periodic orbits converge as
$T\rightarrow\infty$. At this stage it is convenient to note that the Floquet
exponents are the period averages of the “local exponents” Bosetti and Posch
(2014); Ding  _et al._ (2016)

$\lambda_{k}(t)=\frac{\mathbf{w}_{k}(t)^{\top}[\mathbf{f}_{\mathbf{u}}^{\top}(t)+\mathbf{f}_{\mathbf{u}}(t)]\mathbf{w}_{k}(t)}{\|\mathbf{w}_{k}(t)\|^{2}},$
(33)

uniquely defined functions of state space Trevisan and Pancotti (1998)
expressing the local growth rate of tangent perturbations along the invariant
subspaces (here $\|\cdot\|$ indicates the euclidean norm). By the Central
Limit Theorem, the distribution of the $k$-th Floquet exponent across distinct
orbits of similar period $T$, must converge in law to a Dirac delta function
Ott (2002) with standard deviation decaying as $T^{-1/2}$, assuming that the
auto-correlation of time histories of (33) decays sufficiently quickly Ott
(2002). Hence, the bound (32) indicates that the distribution of the
sensitivity of period averages $\mathrm{d}\mathcal{J}/\mathrm{d}\alpha$ will
also converge to a delta function as $T$ increases. In other words, while some
scatter might be observed for short cycles, long periodic orbits will
asymptotically provide the same sensitivity to parameter perturbations.

##  VI Numerical results

To answer the question posed in the introduction, we now turn to numerical
experiments and consider periodic orbits of two well-known chaotic systems.
The first is given by the Lorenz equations Lorenz (1963); Sparrow (2012)

$\begin{cases}\mathrm{d}u_{1}/\mathrm{d}t&=\,\sigma(u_{2}-u_{1}),\\\
\mathrm{d}u_{2}/\mathrm{d}t&=\,\rho u_{1}-u_{2}-u_{1}u_{3},\\\
\mathrm{d}u_{3}/\mathrm{d}t&=\,u_{1}u_{2}-\beta u_{3},\end{cases}$ (34)

where standard parameters $\sigma=10$, $\beta=8/3$ and $\rho=28$ are used
throughout. As in other sensitivity studies on the Lorenz equations (Lea  _et
al._ , 2000; Eyink  _et al._ , 2004; Wang, 2013; Liao, 2016; Lasagna  _et al._
, 2019; Craske, 2019), we consider the sensitivity of the period average of
the observable $J(t)=u_{3}(t)$ with respect to perturbations of $\rho$.
Numerical integration of chaotic trajectories is performed using a classical
fourth-order Runge-Kutta method with $\Delta t=0.005$.

The second system is a finite-dimensional truncation of a spatially extended
system, the Kuramoto-Sivashinky equation

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}+\frac{\partial^{2}u}{\partial x^{2}}+\nu\frac{\partial^{4}u}{\partial
x^{4}}=0,$ (35)

defining the evolution of zero-mean, spatially-periodic fields $u\equiv
u(x,t)$ over the domain $x\in[0,2\pi]$ with dynamics restricted to the
invariant subspace of odd solutions Lan and Cvitanović (2008). Here, we
consider a relatively high diffusivity constant $\nu=(2\pi/L)^{2}$, with
$L=39$, the same value we considered in previous work Lasagna (2018). The
spectral expansion

$\displaystyle u(x,t)=\sum_{k=-n}^{n}iu_{k}(t)\exp(ikx),$ (36)

with $u_{k}=-u_{-k}$, $u_{0}=0$, is truncated at $n=28$, leading to a system
of ODEs approximating solutions of the original partial differential equation.
Some of the numerical results reported in the next sections have been checked
in a statistical sense at finer resolutions, with negligible quantitative
changes. We take the energy density

$J[u(x,t)]=\frac{1}{4\pi}\int_{0}^{2\pi}u^{2}(x,t)\,\mathrm{d}x$ (37)

as the functional of interest and examine the sensitivity of its average with
respect to the diffusivity $\nu$. Numerical integration of chaotic
trajectories is performed using the fourth-order accurate implicit-explicit
method IMEXRKCB4 Cavaglieri and Bewley (2015), with time step $\Delta
t=0.125\nu$.

###  VI.1 Search of long periodic orbits

We use a global Newton-Raphson search algorithm developed in previous work
Lasagna (2018), based on the original method of Ref. [Lan and Cvitanović,
2004] and classical techniques for nonlinear boundary value problems Ascher
_et al._ (1995). Briefly, at iteration $k$, this method solves a Newton-
Raphson update equation to adjust a trial solution composed of a state-space
loop $\mathbf{u}_{k}\in\mathcal{P}^{n}_{T_{k}}$ and a period $T_{k}$. The loop
is not, at least initially, a solution of the equations, i.e. the residual

$\mathbf{r}_{k}(t)=\mathrm{d}{\mathbf{u}}_{k}(t)/\mathrm{d}t-\mathbf{f}(\mathbf{u}_{k}(t),\alpha)\;\in\mathcal{P}_{T_{k}}^{n}$
(38)

is generally different from zero along the loop. The only significant
modification that we have implemented is that the loop derivative operator
$\mathrm{d}/\mathrm{d}t$ is approximated using an eight-order accurate finite-
difference stencil (instead of fourth order), enhancing the overall
accuracy/cost ratio and allowing longer orbits to be found. The same high-
order discretisation is used for the solution of the tangent problem (24).
Initial guesses are obtained from near recurrences of the chaotic flow.

In previous work Lasagna (2018), we attempted to locate exhaustively all short
periodic orbits and examined their sensitivity as a function of the
topological length. In this work, we adopt a different strategy, motivated by
the objective of examining the properties of typical long periodic orbits
found by Newton-Raphson searches. To this end, rather than considering the
topological length, we select a number of arbitrary reference periods, denoted
as $\hat{T}$, and locate up to five thousand periodic orbits with actual
period falling within $\pm 5$% of the reference period. Periodic orbits do
have inherent time scales (the period of the shortest cycle), but for long
reference periods we have observed that the actual period of converged
solutions is uniformly distributed in the $\pm 5$% range. Hence, the reference
periods can be selected arbitrarily and are chosen here such that the maximum
reference period $\hat{T}_{\mathrm{max}}$ is about two orders of magnitude
larger than the period of the shortest admissible cycle, as indicated in table
1. This range is sufficiently wide to reveal the asymptotic convergence of
properties of long periodic orbits as the period increases.

Search results are reported in figure 1.

Figure 1: Number of unique periodic orbits found as a function of the
reference period $\hat{T}$, for the Lorenz equations, panel (a), and KS
system, panel (b). The red line is at 5000. For the KS system, the reference
period is scaled by the diffusivity.

For short reference periods, we locate as many orbits as it is feasible and
the number of periodic orbits found grows exponentially with the period. The
number of periodic orbits found quickly saturates the set threshold of five
thousand orbits. We stress the fact that our focus is not to provide a
description of the statistical distribution of the complete set of periodic
orbits of high topological lengths, nor to use such quantities to approximate
the measure using cycle averaging theory. Such calculations would be biased by
the sampling. Rather, we aim to develop an understanding of what to expect
from long orbits that might be typically found numerically by converging near
recurrence events. Hence, the threshold is chosen so that statistics over the
ensemble are sufficiently robust.

We report in figure 2 the shortest and longest periodic orbits found for the
Lorenz equations, panel (a), and KS system, panel (b). While the short cycles
are topologically simple, longer orbits wind around the attractor in a
complicated fashion and are thus indistinguishable to the eye from long
chaotic trajectories.

Figure 2: Projections of the shortest (solid red) and longest (thin black) periodic orbits found in this study for the Lorenz, panel (a), and KS equations, panel (b). Table 1: Reference temporal grid spacing for the finite-difference approximation of the derivative involved in the search of periodic orbits, period of the shortest admissible orbit and ratio between the largest reference period and $T_{\mathrm{min}}$. | reference $\Delta t$ | $T_{\mathrm{min}}$ | $\hat{T}_{\mathrm{max}}/T_{\mathrm{min}}$
---|---|---|---
Lorenz | 0.01 | 1.5586 |  $\sim$ 645
KS | 0.125$\nu$ | 24.9080$\nu$ |  $\sim$ 201
Figure 3: Convergence of the norm of the residual (38) during the Newton-
Raphson search of long periodic solutions of the Lorenz equation, panel (a),
and KS equations, panel (b), for one hundred test searches with reference
period $\hat{T}=1000$, and $\hat{T}=5000\nu$, respectively. Only converged
searches are reported. The horizontal red line indicates the stopping
tolerance. Panels (c) and (d): success rate of the search of periodic orbits
estimated from the same test searches, as a function of the reference period.

We build an inventory of unique solutions by ensuring that all periodic orbits
found have period average (16) differing to at least eight decimal places.
This is few orders of magnitude larger than the accuracy at which this
quantity is determined in the search, with the temporal discretisation
settings reported in table 1. Because of the high number of orbits existing at
such high topological lengths, very few duplicates have been found. For both
systems, we have observed that relatively few Newton-Raphson iterations are
required and that the convergence history is independent of the period $T$.
However, the success rate appears to decline slightly for longer periods,
arguably as a result of the increasing condition number of the Newton-Raphson
update problem. This is illustrated in figure 3, where we show the history of
the norm of the residual (38) for about a hundred long period searches for the
Lorenz and KS equations, in panels (a) and (b), respectively. Panels (c) and
(d) shows how the success rate, estimated from these test searches, varies
with the period. We have observed these trends to be independent of the
temporal discretisation and, for the KS equations, of the spatial resolution.
This gives us enough confidence that these orbits are numerically reliable
approximations of exact solutions of the equations and not an artefact of the
search method.

In the next sections, we examine properties of these orbits. We first focus on
period averages and Floquet exponents in sections VI.2 and VI.3, respectively.
To characterise how properties vary across the ensemble of orbits at each
reference period, we compute the ensemble mean and standard deviation and
denote these quantities by $\mathrm{mean}[\cdot]$ and $\mathrm{std}[\cdot]$,
respectively. We do not consider these moments as substitutes of cycle
averaging formulae but we use them to characterise in statistical terms the
properties of typical orbits obtained numerically. Probability distributions
of these quantities are also shown, with the caveat that they only represent
typical orbits, and not the complete set of admissible periodic orbits.
Sensitivity of period averages are then finally considered in section VI.4. We
also compare averages, exponents and sensitivities to analogous quantities
computed on chaotic trajectories, to address the original question whether
long periodic orbits can be considered as accurate proxies for the chaotic
state.

###  VI.2 Statistics of time averages

The statistics of period averages over typical long periodic orbits are
compared to statistics of time averages of chaotic trajectories of same
reference length $\hat{T}$ in figure 4. Here and in subsequent figures, error
bars denote plus/minus three times the standard error Rao (1973).

Figure 4: Ensemble mean, panels (a, b), and standard deviation, panels (c, d),
of the time averaged observable as a function of the reference period
$\hat{T}$, for chaotic trajectories (filled diamonds) and periodic orbits
(open circles). In Left panels, Lorenz equations, right panels, KS system.

The notation $\mathcal{J}^{\hat{T}}$ emphasizes the dependence of averages on
the reference period. Panels (a) and (b) show the mean period average of
typical periodic orbits, for the Lorenz and KS systems, compared to the long
time average from chaotic simulations. For periodic orbits, we observe a small
positive/negative distortion at low periods for the KS/Lorenz equations, in
agreement with previous work Zaks and Goldobin (2010). This distortion decays
as $\hat{T}$ increases. More importantly, in statistical terms, the period
average of typical long periodic orbits spanning increasingly larger fractions
of the attractor appears to converge to the long-time average of chaotic
trajectories.

The standard deviation of the finite-time average over chaotic trajectories,
panels (c) and (d), decays asymptotically as $\hat{T}^{-1/2}$. This is the
trend predicted for an ensemble of averages by the Central Limit Theorem Ott
(2002) if correlations decay sufficiently fast. In fact, both systems
considered here display ‘typical’ chaos (in the terminology of Ref. [Prasad
and Ramaswamy, 1999]), with correlations dying out exponentially.
Asymptotically, the standard deviation of the period average decays in the
same manner, i.e. longer periodic orbits provide more accurate descriptions of
the long-time mean. On the other hand, for short periods, the standard
deviation of chaotic trajectories decays at a faster rate for both systems
considered. We argue that this is an effect of correlations affecting the
asymptotic behaviour. For short periodic orbits, increasing $\hat{T}$ does not
necessarily result in a lower variance, in agreement with previous
observations Saiki and Yamada (2009); Zaks and Goldobin (2010). Overall, the
data in figure 4 suggests that typical long periodic orbits found in
computations have, in statistical terms, similar properties to those of
chaotic trajectories.

Figure 5: Probability distributions of time averages over periodic
trajectories, panels (a) and (b), and chaotic trajectories, panels (c) and
(d), for two different reference time spans $\hat{T}$. Left panels, Lorenz
equations, right panels, KS system.

Probability distributions of time averages over periodic and chaotic
trajectories are reported in figure 5 for two different $\hat{T}$. Data for
periodic and chaotic trajectories is reported in panels (a-b) and (c-d),
respectively. Data for the Lorenz equations and the KS system is reported in
the left and right panels, respectively. Before computing the distributions,
samples are normalized such as to have zero mean and unit variance, and the
empirical distributions are compared to the normal distribution, indicated in
grey in the figure. The distributions are skewed to positive values for short
reference time spans. For the Lorenz equations this was observed in previous
work Saiki and Yamada (2009). Eventually, however, the distributions of
averages over chaotic and periodic trajectories appear to collapse to the
normal distribution, although less pronouncedly for periodic orbits. Similar
to the trends of the standard deviations in figure 4, this is the behaviour
dictated by the Central Limit Theorem Ott (2002), for which the probability
distribution of time averaged quantities can asymptotically be approximated
near its peak and within few standard deviations by a Gaussian probability
distribution. These results support the observations reported in Ref. Zaks and
Goldobin (2010), in which the distributions of period averages converge to
delta functions as $\hat{T}\rightarrow\infty$.

###  VI.3 Statistics of Floquet exponents

We now compare statistics of Floquet exponents of typical periodic orbits with
finite-time Lyapunov exponents (FTLE) calculated over chaotic trajectories
with same reference period. We calculate the FTLEs using classical methods
Benettin  _et al._ (1980); Shimada and Nagashima (1979), involving propagating
a set of vectors in tangent space and occasionally performing a re-
orthogonalization using the Gram-Schmid procedure to counter the inevitable
alignment to the most unstable subspace.

Figure 6: The real part of the fifteen leading Floquet exponents of a long and
a short orbit of the KS system. The inset shows the eight leading exponents.
The ninth exponent has a much lower value, at around -82. Figure 7: Evolution
of mean (top panels) and standard deviation (bottom panels) of few selected
Floquet (open circles) exponent and FTLEs (filled diamonds) as a function of
the reference period for the Lorenz equations, panels (a) to (d), and for KS
system, panels (e) to (l).

Computing the spectrum of Floquet exponents is, however, a notoriously
challenging problem Lust (2001); Ding and Cvitanović (2016), even for orbits
of moderate period. The exponential growth of the entries of the monodromy
matrix (6) makes its eigenvalue decomposition inaccurate in finite-precision
arithmetic. In this work, where the focus is on long periodic solutions, we
have used a more robust algorithm recently discussed in Ref. [Ding and
Cvitanović, 2016]. With this algorithm we have been able to compute accurately
Floquet exponents for arbitrarily long orbits, corresponding to multipliers
spanning thousands of orders of magnitude. Our implementation does not make
use of iterative QR-based eigenvalue algorithms, but works in a matrix-free
fashion and only requires an existing time-stepper code for the linearised
equations. The algorithm and a small addition to Ref. Ding and Cvitanović
(2016) are presented in appendix A for completeness.

For the KS system, all periodic orbits found have only one unstable Floquet
eigendirection, and with the present spatial resolution twenty-six contracting
directions. For illustrative purposes, the leading part of the Floquet
spectrum of one long and one short periodic orbit is shown in figure 6. The
first eight exponents fall in the range $(-30,5)$, while the ninth exponent is
sharply more negative and is followed by a long tail of negative exponents.
These correspond to contracting “spurious” modes Ding  _et al._ (2016); Yang
_et al._ (2009), with a value that is closely determined by the linear term of
the governing equations. All orbits in our database have a similar spectrum.

We report in figure 7 data for the two non-trivial exponents of the Lorenz
equations, panels (a) to (d) and the first three non-trivial exponents for the
KS system, panels (e) to (l). The evolution of the mean (top five panels) and
standard deviation (bottom five panels) of selected Floquet exponents (open
circles) and FTLEs (filled diamonds) is reported as a function of the
reference period $\hat{T}$. Similar to figure 4, error bars define plus/minus
three times the standard error. These are shown only for Floquet exponents,
since statistics of the FTLEs are computed over a sufficiently large
collection of independent orbits to make the bars smaller than the symbols in
the figure. As the reference period increases, the average Floquet exponents
of typical orbits converge, within the statistical relevance of our ensemble,
to the corresponding infinite-time Lyapunov exponents. In other words, typical
long periodic orbits found using Newton-Raphson searches have the same
stability properties of long ergodic trajectories. Similar to period averages,
we observe a distortion over short periodic trajectories. The standard
deviation of exponents of periodic and chaotic trajectories decays as
$\hat{T}^{-1/2}$. For short periods, the standard deviation of FTLEs of the
Lorenz equations decreases more rapidly, as $\hat{T}^{-1}$. This is induced by
exponential tails characterizing the distribution of short-time FTLEs, often
observed for intermittent systems Prasad and Ramaswamy (1999).

Figure 8: Probability distributions of the two non-trivial Floquet exponents
and FTLE for the Lorenz equations at $\hat{T}=1000$, panels (a) and (b), and
the first three non-trivial exponents of the KS system at $\hat{T}=5000$,
panels (c) to (e). The parabolas denote Gaussian distributions with mean and
standard deviation equal to the sample mean and variance of the numerical
data. Figure 9: Panel (a): probability distributions of the gradient
$\mathcal{J}_{\rho}^{\hat{T}}$ for typical long periodic orbits of the Lorenz
equations. Panel (b): $(u_{1},u_{3})$ projection of the orbit
$\mathrm{A}^{14}\mathrm{B}$ with period $T=10.4965...$ and extremal gradient
$\mathcal{J}_{\rho}=0.9988...$, for two values of $\rho$.

Probability distributions of the exponents are reported in figure 8, where we
present data for the longest reference period considered to illustrate the
asymptotic behaviour. Normal distributions with mean and variance equal to the
sample mean and variance of the numerical data are reported in grey. Since
Floquet exponents and FTLEs are averages of local quantities (33), their
distributions follow the same trend as that of the quantity
$\mathcal{J}^{\hat{T}}$. Hence, for asymptotically long periods, the
distributions collapse to the same normal behaviour, at least within
plus/minus five standard deviations shown in the figure.

In summary, the distributions of the Floquet exponents of typical periodic
orbits found in our computations localise around the Lyapunov exponents of the
chaotic flow. Based on the discussion of section V and the bound (32), this
suggests that the distribution of the sensitivity of period averages of
typical long orbits will also localise around an asymptotic value. This
localisation is examined in the next section.

###  VI.4 Statistics of sensitivities

The probability distribution of gradients $\mathcal{J}_{\rho}^{\hat{T}}$ from
periodic orbits of the Lorenz equations is reported in figure 9-(a), for three
reference periods. For short periods, the probability distributions are not
normal, but have a heavier left tail that decays fast enough for the first and
second central moments to be finite. Extremal periodic orbits with low
sensitivity in the left tails feature close passes to the unstable equilibrium
at the origin Eckhardt and Ott (1994). Continuation in $\rho$ of the extremal
orbit found for $\hat{T}=10$ (with symbol sequence $\mathrm{A}^{14}\mathrm{B}$
in the notation of Ref. [Viswanath, 2003]), reported in panel (b), shows that
points near the origin of state space move towards the origin, thus causing a
lower sensitivity of the $u_{3}$ variable. For the longer orbits found in our
computations from near recurrences, we observe that the fraction of the period
spent in the neighbourhood of the origin diminishes in relative terms and
approaches that of long chaotic trajectories. As a result, the left tail of
the distributions in figure 9-(a) shows a progressively faster decay as
$\hat{T}$ increases and the distribution ultimately converges to a Gaussian
law (denoted with a dashed line), localised around
$\mathcal{J}_{\rho}^{\hat{T}}\simeq 1.017$.

Figure 10: Panel (a): probability distributions of the gradient
$\mathcal{J}_{\nu}^{\hat{T}}$ for periodic orbits of the KS system. Panel (b):
distributions of the deviation from
$\mathrm{mean}[\mathcal{J}_{\rho}^{\hat{T}}]$ for $\hat{T}=500$. The grey
parabola is a Gaussian fit to the distribution for $\hat{T}=5000\nu$.

Probability distributions of the gradient $\mathcal{J}^{\hat{T}}_{\nu}$ for
typical periodic orbits of the KS system are reported in panel (a) of figure
10. Near the peak, the distributions can be approximated reasonably well by a
Gaussian law. However, much higher/lower sensitivities are observed for a few
orbits, resulting in a significant departure from normality and heavy tails.
To characterize these tails more precisely, we find fifteen thousand more
periodic orbits for $\hat{T}=500$ and show in panel (b) the probability
distribution of the deviation from the mean of the distribution in panel (a),
for this reference period. The tails are well described by a power-law
distribution of the form $p(x)=x^{-n}$ with exponent $n=3$. This structure is
an inherent feature of the problem and not a numerical artefact depending, for
instance, on the resolution.

As illustrated in section V, large sensitivities can be directly associated to
bifurcations. This is illustrated in figure 11, where we report the
continuation analysis of an orbit at $\hat{T}=500$ with large gradient
$\mathcal{J}_{\nu}\simeq 1572.27$. The average energy density
$\mathcal{J}(\nu)$ and its gradient $\mathcal{J}_{\nu}(\nu)$ are reported as a
function of the bifurcation parameter $\nu$ in panels (a) and (b),
respectively.

Figure 11: Panel (a): period average of an orbit with $\hat{T}=500$ with large
positive gradient $\mathcal{J}_{\nu}\simeq 1572.27$, continued over the
parameter $\nu$. Panel (b): sensitivity of the period average for the same
orbit. Open circles denote data from the continuation. The red curve is the
model (40) fitted to the data. The smaller inset shows the same quantity on
shifted coordinates, in a bi-logarithmic plot.

Near the bifurcation point, denoted by $\nu_{b}$, the period average is well
described by the functional form

$\mathcal{J}(\nu)\simeq
C_{0}+C_{1}\sqrt{\nu-\nu_{b}}+C_{2}(\nu-\nu_{b}),\quad\nu>\nu_{b},$ (39)

where the square root term is typical in normal forms of bifurcations for
periodic orbits Guckenheimer and Holmes . The gradient is then

$\mathcal{J}_{\nu}(\nu)\simeq\frac{C_{1}}{2\sqrt{\nu-\nu_{b}}}+C_{2}$ (40)

and approaches infinity at $\nu_{b}$ as $1/\sqrt{\nu}$. Fitting the data in
panel $(b)$ to the model (40) shows that the constant $C_{2}$, the gradient
measured sufficiently far away from the bifurcation point, is about $-246.65$,
in line with the high probability region of the distributions in figure
10-(a).

The functional form (39) is sufficient to explain the structure of the tails
in figure 10. To this end, assume that periodic orbits appear in bifurcations
at critical values $\nu_{b}$ as $\nu$ is increased, as expressed by (39).
Focusing at a given $\hat{T}$, assume also that the number of periodic orbits
is large, so that $\mathcal{J}(\nu)$ can be thought of as a random variable,
with the coefficients $C_{0},C_{1},C_{2}$ and the bifurcation point $\nu_{b}$
being random variables with values differing from orbit to orbit. The gradient
$\mathcal{J}_{\nu}(\nu)$ is then also a random variable that can take
arbitrarily large values if $\nu-\nu_{b}$ is small. Now, the probability that
the gradient $\mathcal{J}_{\nu}(\nu)$ is less than some large positive
constant $x$ can be expressed by introducing the cumulative distribution
function $P_{\mathcal{J}_{\nu}}(x)$, defined as

$\displaystyle P_{\mathcal{J}_{\nu}}(x)=\;$
$\displaystyle\mathrm{prob}\left[\mathcal{J}_{\nu}(\nu)<x\right]$
$\displaystyle=\;$ $\displaystyle
1-\mathrm{prob}\left[\mathcal{J}_{\nu}(\nu)>x\right]$ $\displaystyle=\;$
$\displaystyle 1-\mathrm{prob}\left[\nu-\nu_{b}<(C_{1}/2x)^{2}\right],$ (41)

where we have used the definition (40) and neglected $C_{2}$, since $x\gg 1$,
to develop the algebra in the last step. The probability in the third line can
be equivalently interpreted as the probability that bifurcation points are
closer to the reference value than a distance $(C_{1}/2x)^{2}$. Assuming the
points $\nu_{b}$ not to be preferentially distributed on the real line near
$\nu$, this probability is then $cx^{-2}$ for some constant $c$. In other
words, the larger $x$, the less likely is that a periodic orbit bifurcates
near $\nu$. Hence, the cumulative distribution of $P_{\mathcal{J}_{\nu}}(x)$
must, asymptotically for large $x$, obey

$P_{\mathcal{J}_{\nu}}(x)=1-cx^{-2}.$ (42)

The probability distribution of the gradient $\mathcal{J}_{\nu}$ can then be
obtained by differentiating the cumulative distribution with respect to its
argument, leading to the power-law

$p(x)=cx^{-3},\quad x\gg 1,$ (43)

the behaviour observed in figure 9. More generally, sampling functions that
have poles of the form $(\nu-\nu_{b})^{\gamma}$, produces probability
distributions with power-law tails of the form $p(x)\simeq x^{-n}$, $x\gg 1$,
with exponent $n=(\gamma-1)/\gamma$, leading to $n=3$ for the present case
with $\gamma=-1/2$ of equation (40).

Figure 12: Dependence of the average, panels (a) and (b), and standard
deviation, panels (c) and (d), of the gradient $\mathcal{J}^{\hat{T}}_{\rho}$,
for the Lorenz equations (left panels), and $\mathcal{J}^{\hat{T}}_{\nu}$ for
the KS system (right panel), on the reference period $\hat{T}$. For the KS
equations, the median and interquartile range are used, denoted as
$\mathrm{median}[\cdot]$ and $\mathrm{iqr}[\cdot]$. Error bars for the KS
system are estimated using a bootstrapping technique.

We now come to the central result of this paper and examine the sensitivity of
typical periodic orbits as a function of the reference period. For
distributions with power-law tails of the form $p(x)\simeq x^{-n}$, central
moments of order $m$ are undefined for $m\geq n-1$. For the KS system, with
$n=3$, while the mean sensitivity across periodic orbits is defined (although
convergence is weak), the standard deviation is not. Hence, for this system,
we use the median and interquartile range, as measures for the localisation
and variability of gradients, respectively Results are shown in figure 12.
Panels (a) and (b) show the mean and median gradient for the Lorenz and KS
systems, respectively. Panels (c) and (d) show the standard deviation and
interquartile range of the gradient. We observe that, as
$\hat{T}\rightarrow\infty$ the mean/median converges to a value that is
approximately $\mathcal{J}_{\rho}$ = 1.017 and $\mathcal{J}_{\nu}=-155$, for
the Lorenz and KS systems, respectively. Given the bound (32), the convergence
of sensitivities is consistent with the convergence of the Floquet exponents
of figure 7. However, for short cycles, the sensitivity of periodic orbits can
be, in average terms over the inventory of available orbits, remarkably
different to that of long cycles. This behaviour is more pronounced for the KS
system. The standard deviation and interquartile range of the sensitivity
follow the same asymptotic behaviour of the period averages and decay
asymptotically as ${\hat{T}}^{-1/2}$, indicating that the corresponding
probability distributions localise around the averages of panels (a) and (b).
An important remark is that sensitivity computations using periodic orbits do
not suffer from shadowing errors displayed by shadowing methods applied to
chaotic trajectories, e.g. the Least-Squares Shadowing Wang (2014) and
Periodic Shadowing algorithms Lasagna  _et al._ (2019). For such algorithms,
convergence proofs have been offered that suggest that the standard deviation
of sensitivity calculations on hyperbolic systems should first decay as
${\hat{T}}^{-1}$ as a result of the approximations of the exact shadowing
direction involved in these algorithms.

Overall, these results show that the sensitivity computed from typical long
periodic orbits found in computations converges to a well defined value as the
period increases. This is not to say that all long orbits can be considered
good proxies, as the tails of figure 10 demonstrate, but rather that
sensitivities computed from longer orbits found in computation are likely to
be closer to the asymptotic value. This asymptotic value of the sensitivity is
now compared with the response of long-time averages to finite-amplitude
parameter perturbations using long chaotic simulations. Carefully conducted
numerical approximations of the gradient using a finite-difference formula
(see Ref. [Lasagna  _et al._ , 2019] for details) show that the response of
the average of $u_{3}$ to perturbation of $\rho$ in the Lorenz equations is
approximately $\mathcal{J}_{\nu}\simeq 1.002$, well below the asymptotic value
from long periodic orbits. We remark that this difference is likely not a bias
arising from using using periodic orbits. In fact, the same difference has
been previously observed using other shadowing algorithms applied to chaotic
trajectories Lasagna  _et al._ (2019). Numerical evidence has been provided
Lucarini (2009) suggesting that the Lorenz equations have a linear response to
perturbations of the parameter $\rho$, despite not being hyperbolic (it is a
singularly hyperbolic system in the terminology of Ref. [Morales  _et al._ ,
1999]). For this system, it has also been speculated Reick (2002) that some
observables might vary continuously with parameters and that the bifurcating
orbits have very long period and their effect of the invariant measure is
negligible. How to reconcile the existence of a linear response with the
difference we observe between the prediction of shadowing methods and the
actual response of the system is a question that deserves further analysis.

For the KS system, we show in figure 13 how the long-time averaged energy
density varies with the diffusivity. The data points are obtained by first
computing averages over tens of thousands of long independent segments of
length $T=5000\nu$, and then reporting the median value, which is more robust
to outliers arising from initial conditions leading to a non-chaotic state.
Using a bootstrapping technique we have also computed the standard error on
the median, which is typically smaller than the symbol size in the figure, and
it is thus not shown. We use this approach, instead of reporting the time
average of one long chaotic computation, as it provides a measure of the
accuracy of long-time average estimate. Panel (b) and (c) focus near the
reference diffusivity $\nu=(2\pi/39)^{2}$ in the area spanned by the vertical
lines in panels (a) and (b), respectively.

Figure 13: Long-time mean of the energy density as a function of $\nu$. Data
points denote the median time average across thousand of simulations from
different initial conditions, with averaging time $T=5000\nu$. The dashed line
represents the slope predicted by periodic orbits with reference period
$\hat{T}=5000\nu$. Panels (b) and (c) focus on the area between the two red
vertical lines in panels (a) and (b), respectively.

The dashed line represents the asymptotic gradient $\mathcal{J}_{\nu}$ from
periodic orbits. The system clearly lacks a linear response, in the sense that
the limit

$\lim_{\delta\nu\rightarrow
0}\frac{\mathcal{J}^{\infty}(\nu+\delta\nu)-\mathcal{J}^{\infty}(\nu)}{\delta\nu},$
(44)

is not defined, as the response of the system is not proportional to the
perturbation in the parameter Reick (2002); Farmer (1985); Ershov (1993), at
any scale. As the distributions of sensitivities in figure 10 suggest, some
orbits are always infinitesimally close to bifurcation and small parameter
perturbations might induce abrupt changes in the structure of the attractor,
making chaotic averages non-differentiable. In such conditions, the meaning of
gradients obtained from linear methods, either on periodic or chaotic
trajectories, is unclear.

##  VII Conclusions

In this paper we set out to address the question whether typical long periodic
orbits found numerically may be used as accurate proxies for the sensitivity
of the chaotic state to parameter perturbations. Our motivation to address
this question arises from well known challenges in locating periodic orbits in
fluid systems governed by the Navier-Stokes equations, and the consequential
possibility that predictions of cycle averaging formulae using an incomplete
set might be inaccurate. If the answer to the above question is affirmative,
an heuristic strategy would be to spend available computational resources to
locate one or few orbits, of sufficiently long period. Accurate sensitivity
information from these orbits may then facilitate control and optimisation
tasks.

Here, we have considered long periodic orbits of two low-dimensional chaotic
systems, the Lorenz equations at standard parameters and a minimal-domain
Kuramoto-Sivashinky system with dynamics restricted to the anti-symmetric
subspace. We have built an inventory of thousands of periodic orbits with
period up to two orders of magnitude larger than the shortest admissible
cycle. This approach was not guided by the idea of obtaining an exhaustive
hierarchy of cycles and a complete understanding of their properties. Rather,
we aimed to examine in statistical terms the properties of typical long
orbits. This analysis is naturally biased by the search process, but the bias
leans precisely towards the direction required to answer our original
question, i.e. it favours orbits that can be found in practical computations.

One conclusion from this study is that period averages of long orbits appear
to converge to the long time average of chaotic trajectories. Floquet
exponents, being the period averages of the local rate of growth of
infinitesimal perturbations, also exhibit the same behaviour. Hence, Floquet
exponents of long orbits converge to the Lyapunov exponents calculated using
standard methods. These results are interesting, but perhaps not so important
from an operational point of view, since these quantities can be calculated
directly from chaotic trajectories at a lower cost. The important result is
that the sensitivity of period averages of typical periodic orbits also
converges to a defined value as the period increases. This result is
consistent with the convergence of the Floquet exponents, based on the
relation between stability and sensitivity of orbit, established in this
paper. Interestingly, the probability distribution of sensitivities from
typical orbits can display power-law tails of the form $p(x)\simeq x^{-3}$.
This result can be explained by using a statistical argument and classical
normal forms of bifurcations of periodic orbits. In practice, this suggests
that sensitivity information from orbits found numerically may occasionally be
quite inaccurate.

Some open challenges are now considered. First, as also observed elsewhere
with other shadowing algorithms Blonigan and Wang (2018); Lasagna  _et al._
(2019) the sensitivity of most typical long periodic orbits is not necessarily
consistent with the response of the system to finite-amplitude parameter
perturbations. In absence of differentiability, the meaning and value of
sensitivities from periodic orbits (or other shadowing algorithms) is unclear
and deserves further analysis. The conjecture is that for high-dimensional
systems, where statistics behave as if a linear response existed Albers  _et
al._ (2006); Ruelle (1999), the “thermodynamic limit” of Ruelle Ruelle (1999),
may be invoked and a better consistency between sensitivities from shadowing
methods and the response of the system might be observed Blonigan and Wang
(2018); Lasagna  _et al._ (2019). However, evidence in support of this
conjecture is currently lacking.

A second challenge is that the applicability of the ideas discussed in this
paper to fluid systems governed by the Navier-Stokes remains unclear (see Ref.
Chandler and Kerswell (2013)). One major issue is that the increase of system
dimension inevitably implies a decrease of good near recurrence events. As
advocated in Ref. van Veen  _et al._ (2019) more robust search methods are
required. A third potential issue is that the size of the linear systems
arising in the Newton-Raphson search iterations grows linearly with the period
$T$, regardless of the numerical method utilized, i.e. either for global
search methods Lan and Cvitanović (2004); Lasagna (2018) or with multiple-
shooting techniques Ascher  _et al._ (1995); Cvitanović  _et al._ (2010). The
condition number of these problems grows with $T$, introducing errors in the
Newton-Raphson corrections that might eventually prevent convergence. Hence,
finding long periodic orbits might, eventually, prove too challenging.

##  Appendix A Floquet exponents of long periodic orbits

The algorithm used to compute Floquet exponents of long periodic orbits
follows closely the approach originally introduced in Ref. [Ding and
Cvitanović, 2016]. For completeness, we describe in this appendix this
algorithm and outline our novel contribution.

The algorithm exploits two fundamental facts. The first is that the Jacobian
matrices (3) obey the multiplicative property

$\mathbf{M}(t_{2},t_{0})=\mathbf{M}(t_{2},t_{1})\mathbf{M}(t_{1},t_{0}),$ (45)

for any times $t_{2}\geq t_{1}\geq t_{0}$. Hence, the monodromy matrix
$\mathbf{M}(T,0)$ can be equivalently expressed as the product of $M$ short-
time Jacobian matrices $\mathbf{M}_{i}=\mathbf{M}(t_{i},t_{i-1})$,
$i=1,\ldots,M$, as

$\displaystyle\mathbf{M}(T,0)=\mathbf{M}_{M}\mathbf{M}_{M-1}\ldots\mathbf{M}_{1},$
(46)

for a partitioning of the interval $[0,T]$ into $M$ sub-intervals specified by
times $0\equiv t_{0}>t_{1}>\ldots>t_{M-1}>t_{M}\equiv T$. Note that the
Jacobian matrices $\mathbf{M}_{i}$ obey the cyclic property
$\mathbf{M}_{M+1}=\mathbf{M}_{1}$.

The second fact is that a well-conditioned eigenvalue revealing decomposition
exists for products of matrices such as (46). This is the periodic real Schur
decomposition Bojanczyk  _et al._ (1992); Hench and Laub (1994), initially
introduced in the context of Floquet analysis in Ref. [Lust, 2001] for the
computation of the multipliers and more recently extended Ding and Cvitanović
(2016); Ding  _et al._ (2016) to compute the eigenfunctions. This
decomposition consists in factorizing the short-time Jacobian matrices using a
set of orthogonal matrices $\mathbf{Q}_{i}$, $i=1,\ldots,M$, satisfying the
cyclic property $\mathbf{Q}_{0}=\mathbf{Q}_{M}$, as

$\mathbf{M}_{i}=\mathbf{Q}_{i}\mathbf{R}_{i}\mathbf{Q}_{i-1}^{\top}$ (47)

where the factors $\mathbf{R}_{i},i=1,\ldots,M-1$ are upper triangular
matrices and $\mathbf{R}_{M}$ is in real Schur form, a block upper-triangular
matrix with either $1\times 1$ and $2\times 2$ blocks on the diagonal, in case
the monodromy matrix possesses pairs of complex conjugate multipliers.

Using these two facts, the monodromy matrix can be expressed in real Schur
form as

$\mathbf{M}(T,0)=\mathbf{Q}_{0}\mathbf{R}_{M}...\mathbf{R}_{2}\mathbf{R}_{1}\mathbf{Q}_{0}^{\top}.$
(48)

The product $\mathbf{R}_{M}...\mathbf{R}_{2}\mathbf{R}_{1}$ and the monodromy
matrix are unitarily similar and thus share the same spectrum of eigenvalues.
However, because of the structure of the factors $\mathbf{R}_{i}$, obtaining
the spectrum is a straightforward computation, since the spectrum of a block
triangular matrix is the union of the spectra of the blocks. The structure of
the block upper triangular factor $\mathbf{R}_{M}$ determines whether
exponents are real or form complex conjugate pairs (see Ref. [Golub and Van
Loan, 2012], Th. 7.4.1). For a $1\times 1$ block at location $(i,i)$, a real
Floquet exponent can be obtained as

$\lambda_{i}=\log(\mu_{i})/T=\frac{1}{T}\log\prod_{j=1}^{M}[\mathbf{R}_{j}]_{ii}=\frac{1}{T}\sum_{j=1}^{M}\log[\mathbf{R}_{j}]_{ii}$
(49)

Computing the sum of the logarithms is recommended, as multiplication can
quickly over/underflow before the logarithm is taken. For a $2\times 2$ block,
a pair of complex conjugate exponents can be obtained with a bit more work by
recursively multiplying all $2\times 2$ blocks of the factors $\mathbf{R}_{j}$
at location $(i,i+1)$, and accumulating the sum of the logarithms of scaling
factors required to set the largest element in the partial products to have
unitary magnitude. Overall, this algorithm only operates on well-conditioned
short-time Jacobian matrices, instead of forming the monodromy matrix, and it
is numerically robust.

The numerical algorithm required to obtain the factors $\mathbf{R}_{i}$ and
$\mathbf{Q}_{i}$ in equation (47) from the short-time jacobian matrices is
based on classical QR-based eigenvalue algorithms Kressner (2006); Golub and
Van Loan (2012). Developing a robust implementation is a lengthy and delicate
task. In this paper, we have adopted a matrix-free algorithm introduced in
Ref. [Ding and Cvitanović, 2016], which is a specialization to periodic orbits
of classical methods to compute Lyapunov exponents Benettin  _et al._ (1980).
The approach only requires the action of these matrices on a set of tangent
vectors and is thus computationally more efficient when only a handful of the
leading Floquet exponents is required. The algorithm is simpler to implement,
does not require advanced linear algebra technique and only requires minor
modifications to an existing time-stepper code for the linearised equations.

The algorithm is iterative and we denote quantities at iteration $k$, with a
superscript $(k)$. First, a set of $m$ linearly independent tangent vectors is
defined, where $m$ is the number of desired exponents. For notational
convenience, we arrange them as the columns of matrix
$\hat{\mathbf{Q}}_{0}^{k}\in\mathbb{R}^{n,m}$, where we use a hat to denote a
matrix with reduced dimensions. The period is divided into $M$ sub-intervals,
which need not have the same size. In each sub-interval, the columns of
$\hat{\mathbf{Q}}_{0}^{k}$ are: 1) propagated forward in time using a
linearised time-stepping solver and 2) re-orthogonalized in place using a
Gram-Schmidt procedure. These two steps are formally equivalent to computing

$\mathbf{M}_{i}\hat{\mathbf{Q}}_{i-1}^{k}=\hat{\mathbf{Q}}_{i}^{k}\hat{\mathbf{R}}_{i}^{k}\quad
i=1,\ldots,M,$ (50)

which is akin to (47). Note that forming $\mathbf{M}_{i}$ is not necessary,
only its action on the columns of $\hat{\mathbf{Q}}_{i-1}^{k}$ is required,
making the approach suitable for PDE problems. The triangular factors
$\hat{\mathbf{R}}_{i}^{k}\in\mathbb{R}^{m\times m}$ from the orthogonalization
are stored for later processing. The time between subsequent re-
orthogonalizations depends on the expanding/contracting characteristics of the
tangent space and should be chosen such that the columns of
$\hat{\mathbf{Q}}_{i}^{k}$ remain numerically linearly independent.

After the last sub-interval, the iterations are restarted by setting
$\hat{\mathbf{Q}}_{0}^{k+1}=\hat{\mathbf{Q}}_{M}^{k}$ and after $k$ iterations
the monodromy matrix is formally equivalent to

$\mathbf{M}(T,0)=\hat{\mathbf{Q}}_{0}^{k+1}\hat{\mathbf{R}}_{M}^{k}\ldots\hat{\mathbf{R}}_{2}^{k}\hat{\mathbf{R}}_{1}^{k}\hat{\mathbf{Q}}_{0}^{k\top}.$
(51)

If the first $m$ Floquet multipliers are all real and distinct, the columns of
$\hat{\mathbf{Q}}_{0}^{k+1}$ converge geometrically to a basis for the
subspace spanned by the leading $m$ Floquet eigenvectors. Hence, the
difference $\|\hat{\mathbf{Q}}_{0}^{k+1}-\hat{\mathbf{Q}}_{0}^{k}\|$ converges
to zero in any norm and (51) is the real Schur form of the monodromy matrix,
as in equation (48). In fact, this iteration procedure is a form of subspace
iteration Saad (2011) (also known as orthogonal iteration Golub and Van Loan
(2012), and referred to as simultaneous iteration in Ref. [Ding and
Cvitanović, 2016]). The leading $m$ Floquet exponents can then be obtained
from the diagonals of the factors $\hat{\mathbf{R}}_{i}$, as discussed.

Figure 14: The base ten logarithm of the entries of the rotation matrix
$\hat{\mathbf{D}}_{k}$ at iterations $k=1,3,10$ and 33, panels (a-d).
Convergence of error on the estimate of a few selected Floquet exponents,
panel (e).

A simple adjustment of this approach can be introduced when some of the
multipliers form complex conjugate pairs. The iterations still converge, in
the sense that the subspace spanned by $\hat{\mathbf{Q}}_{0}^{k}$ converges,
but only the columns associated to real exponents converge individually Ding
and Cvitanović (2016). The subspace spanned by a pair of columns of
$\hat{\mathbf{Q}}_{0}^{k+1}$ corresponding by the space spanned by the Floquet
eigenvector associated to complex conjugate multipliers also converges, but at
every iteration the two columns are rotated by an angle in the subspace they
span. Hence, we introduce a rotation matrix $\hat{\mathbf{D}}^{k}$ such that

$\hat{\mathbf{Q}}_{0}^{k+1}=\hat{\mathbf{Q}}_{0}^{k}\hat{\mathbf{D}}^{k}.$
(52)

For large $k$, this rotation converges to a matrix that has the structure of
the product of Givens rotation matrices, each rotating one pair of columns of
$\hat{\mathbf{Q}}_{0}^{k}$ to the corresponding pair of
$\hat{\mathbf{Q}}_{0}^{k+1}$. With this modification, the product
$\hat{\mathbf{D}}^{k}\hat{\mathbf{R}}^{k}_{M}$ converges to a block upper
triangular matrix and (51) with (52) is formally equivalent to (48).

To the best of the author’s understanding, a procedure to compute this
rotation matrix was not proposed Ref. [Ding and Cvitanović, 2016], which
focused instead on using iterative QR-based algorithms. The new contribution
of this appendix is a simple strategy to obtain it. The rotation
$\hat{\mathbf{D}}^{k}$ can be found as the solution of the orthogonal
Procrustes problem Golub and Van Loan (2012)

$\underset{\mathbf{D}^{k}}{\mathrm{argmin}}\;\;\|\hat{\mathbf{Q}}_{0}^{k+1}-\hat{\mathbf{Q}}_{0}^{k}\mathbf{D}^{k}\|_{\mathrm{F}}=\mathbf{U}^{k}\mathbf{V}^{k\top},$
(53)

where $\|\cdot\|_{\mathrm{F}}$ is the Frobenius norm and where the two
matrices at the right hand side are obtained from the Singular Value
Decomposition

$\hat{\mathbf{Q}}_{0}^{k\top}\hat{\mathbf{Q}}_{0}^{k+1}=\mathbf{U}^{k}\mathbf{\Sigma}^{k}\mathbf{V}^{k\top}.$
(54)

In our implementation, we compute the rotation $\hat{\mathbf{D}}^{k}$ along
the iterations and use a simple heuristic to detect pairs of complex conjugate
eigenvectors, or inverse hyperbolic directions, when a diagonal entry is close
to $-1$. We then monitor the maximum absolute difference between estimates of
the Floquet exponents and stop the iterations when such difference is lower
than a user-defined tolerance. Panels (a-d) of figure 14 shows the progressive
convergence of the rotation matrix $\hat{\mathbf{D}}_{k}$ for a calculation on
the shortest periodic orbit of the KS system reported in figure 2. Iteration
$k=1,3,10$ and 33 are shown. Except for the fifth and sixth column, the
columns of $\hat{\mathbf{Q}}_{0}^{k+1}$ converge to the columns of
$\hat{\mathbf{Q}}_{0}^{k}$ and all Floquet multipliers are real and positive.
Panel (e) shows the convergence of the error on the estimate of a few Floquet
exponents.

In practice, we have found this method to be quite robust for the systems used
in this paper, where Floquet exponents are typically well separated. We have
observed that the number of iterations required for convergence decreases with
the period, with exponents of the longest periodic orbits of the KS system
requiring only three/four iterations to converge to machine accuracy. This can
be attributed to the faster convergence of the columns of
$\hat{\mathbf{Q}}^{k}_{i}$ to the subspace spanned by the leading Floquet
eigenmodes, as the integration time is proportionally longer, following the
same pattern of convergence of algorithms to compute Lyapunov exponents from
chaotic trajectories Goldhirsch  _et al._ (1987). We have therefore made no
attempt at improving the convergence rate and computational cost by using
shift and deflation techniques that are customarily used in state-of-the-art
implementations of eigenvalue algorithms Golub and Van Loan (2012); Kressner
(2006).

## References

  * Viswanath (2007) D. Viswanath, Journal of Fluid Mechanics 580, 339 (2007).
  * Chandler and Kerswell (2013) G. J. Chandler and R. R. Kerswell, Journal of Fluid Mechanics 722, 554 (2013).
  * Willis  _et al._ (2016) A. P. Willis, K. Y. Short, and P. Cvitanović, Phys. Rev. E 93, 022204 (2016).
  * Budanur  _et al._ (2017) N. B. Budanur, K. Y. Short, M. Farazmand, A. P. Willis, and P. Cvitanović, Journal of Fluid Mechanics 833, 274 (2017).
  * Lasagna (2018) D. Lasagna, SIAM Journal on Applied Dynamical Systems 17, 547 (2018).
  * Hwang and Choi (2008) Y. Hwang and H. Choi, Physics of Fluids (1994-present) 20, 071703 (2008).
  * Giannetti  _et al._ (2010) F. Giannetti, S. Camarri, and P. Luchini, Journal of Fluid Mechanics 651, 319 (2010).
  * Meliga  _et al._ (2016) P. Meliga, E. Boujo, and F. Gallaire, Journal of Fluid Mechanics 800, 327 (2016).
  * Giannetti  _et al._ (2019) F. Giannetti, S. Camarri, and V. Citro, Journal of Fluid Mechanics 864, 45 (2019).
  * Lea  _et al._ (2000) D. J. Lea, M. R. Allen, and T. W. N. Haine, 52, 523 (2000).
  * Luchini and Quadrio (2002) P. Luchini and M. Quadrio, in  _American Society of Mechanical Engineers, Fluids Engineering Division (Publication) FED_ , Universita di Salerno, Salerno, Italy (ASME, 2002) pp. 1381–1385.
  * Bowen (1975) R. Bowen, Journal of differential equations 18, 333 (1975).
  * Ruelle (1999) D. Ruelle, Journal of Statistical Physics 95, 393 (1999).
  * Nusse and Yorke (1988) H. E. Nusse and J. A. Yorke, Commun.Math. Phys. 114, 363 (1988).
  * Wang (2013) Q. Wang, Journal of Computational Physics 235, 1 (2013).
  * Wang  _et al._ (2014) Q. Wang, R. Hu, and P. Blonigan, Journal of Computational Physics 267, 210 (2014).
  * Ni and Wang (2017) A. Ni and Q. Wang, Journal of Computational Physics 347, 56 (2017).
  * Shawki and Papadakis (2018) K. Shawki and G. Papadakis, arXiv (2018), 1810.12222v1 .
  * Eckhardt and Ott (1994) B. Eckhardt and G. Ott, Z. Physik B - Condensed Matter 93, 259 (1994).
  * Cvitanović (1995) P. Cvitanović, Physica D 83, 109 (1995).
  * Artuso  _et al._ (1999) R. Artuso, E. Aurell, and P. Cvitanović, Nonlinearity 3, 325 (1999).
  * Kazantsev (1998) E. Kazantsev, Nonlinear Processes in Geophysics 5, 193 (1998).
  * Kazantsev (2001) E. Kazantsev, Chaos 12, 1989 (2001).
  * Viswanath (2003) D. Viswanath, Nonlinearity 16, 1035 (2003).
  * Lucas and Kerswell (2015) D. Lucas and R. R. Kerswell, Physics of Fluids (1994-present) 27, 045106 (2015).
  * van Veen  _et al._ (2019) L. van Veen, A. Vela-Martín, and G. Kawahara, Phys. Rev. Letters 123, 134502 (2019).
  * Budanur  _et al._ (2015) N. B. Budanur, P. Cvitanović, R. L. Davidchack, and E. Siminos, Phys. Rev. Letters 114, 084102 (2015).
  * Zaks and Goldobin (2010) M. A. Zaks and D. S. Goldobin, Phys. Rev. E 81, 018201 (2010).
  * Lorenz (1963) E. N. Lorenz, J. Atmos. Sci. 20 (1963).
  * Christiansen  _et al._ (1997) F. Christiansen, P. Cvitanović, and V. Putkaradze, Nonlinearity 10, 55 (1997).
  * Lan and Cvitanović (2008) Y. Lan and P. Cvitanović, Phys. Rev. E 78, 026208 (2008).
  * Sánchez  _et al._ (2004) J. Sánchez, M. Net, B. Garcı́a-Archilla, and C. Simó, Journal of Computational Physics 201, 13 (2004).
  * Farkas (1994) M. Farkas,  _Periodic Motions_ , Applied Mathematical Sciences, Vol. 104 (Springer Science & Business Media, New York, NY, 1994).
  * Moore (2006) G. Moore, SIAM J. Numer. Anal. 42, 2522 (2006).
  * (35) J. Guckenheimer and P. Holmes,  _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields_ , Volume 42 of Applied Mathematical Sciences (Spinger-Verlag).
  * Marcotte and Grigoriev (2016) C. D. Marcotte and R. O. Grigoriev, Chaos 26, 093107 (2016).
  * Lasagna  _et al._ (2019) D. Lasagna, A. Sharma, and J. Meyers, Journal of Computational Physics 391, 119 (2019).
  * Viswanath (2001) D. Viswanath, SIAM Rev. 43, 478 (2001).
  * Hale (1969) J. K. Hale,  _Ordinary Differential Equations_ , Vol. 21 (Wiley-Interscience, 1969).
  * Lan and Cvitanović (2004) Y. Lan and P. Cvitanović, Phys. Rev. E 69, 016217 (2004).
  * Eyink  _et al._ (2004) G. L. Eyink, T. W. N. Haine, and D. J. Lea, Nonlinearity 17, 1867 (2004).
  * Ginelli  _et al._ (2007) F. Ginelli, P. Poggi, A. Turchi, H. Chaté, R. Livi, and A. Politi, Phys. Rev. Letters 99, 27 (2007).
  * Wolfe and Samelson (2007) C. L. Wolfe and R. M. Samelson, 59, 355 (2007).
  * Xu and Paul (2018) M. Xu and M. R. Paul, Phys. Rev. E 97, 032216 (2018).
  * Nikitin (2018) N. Nikitin, Journal of Fluid Mechanics 849, 942 (2018).
  * Inubushi  _et al._ (2015) M. Inubushi, S. Takehiro, and M. Yamada, Phys. Rev. E 99, 27 (2015).
  * Luchini and Bottaro (2014) P. Luchini and A. Bottaro, Annu. Rev. Fluid Mech. 46, 493 (2014).
  * Bosetti and Posch (2014) H. Bosetti and H. A. Posch, Commun. Theor. Phys. 62, 451 (2014).
  * Ding  _et al._ (2016) X. Ding, H. Chaté, P. Cvitanović, E. Siminos, and K. A. Takeuchi, Phys. Rev. Letters 117, 024101 (2016).
  * Trevisan and Pancotti (1998) A. Trevisan and F. Pancotti, J. Atmos. Sci. 55, 390 (1998).
  * Ott (2002) E. Ott,  _Chaos in Dynamical Systems_ (Cambridge University Press, 2002).
  * Sparrow (2012) C. Sparrow,  _The Lorenz Equations_ , Bifurcations, Chaos, and Strange Attractors, Vol. 41 (Springer Science & Business Media, New York, NY, 2012).
  * Liao (2016) H. Liao, Journal of Computational Physics 313, 57 (2016).
  * Craske (2019) J. Craske, Chaos, Solitons & Fractals 119, 1 (2019).
  * Cavaglieri and Bewley (2015) D. Cavaglieri and T. Bewley, Journal of Computational Physics 286, 172 (2015).
  * Ascher  _et al._ (1995) U. M. Ascher, R. M. M. Mattheij, and R. D. Russell,  _Numerical solution of boundary value problems for ordinary differential equations_ , edited by SIAM (Society for Industrial and Applied Mathematics, 1995).
  * Rao (1973) C. R. Rao,  _Linear Statistical Inference and its Applications_ , 2nd ed. (John Wiley & Sons, New York, 1973).
  * Prasad and Ramaswamy (1999) A. Prasad and R. Ramaswamy, Phys. Rev. E 60, 2761 (1999).
  * Saiki and Yamada (2009) Y. Saiki and M. Yamada, Phys. Rev. E 79, 015201 (2009).
  * Benettin  _et al._ (1980) G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Meccanica 15, 9 (1980).
  * Shimada and Nagashima (1979) I. Shimada and T. Nagashima, Progress of theoretical physics 61, 1605 (1979).
  * Lust (2001) K. Lust, International Journal of Bifurcation and Chaos in Applied Sciences and Engineering 11, 2389 (2001).
  * Ding and Cvitanović (2016) X. Ding and P. Cvitanović, SIAM Journal on Applied Dynamical Systems 15, 1434 (2016).
  * Yang  _et al._ (2009) H.-L. Yang, K. A. Takeuchi, F. Ginelli, H. Chaté, and G. Radons, Phys. Rev. Letters 102, 1423 (2009).
  * Wang (2014) Q. Wang, SIAM J. Numer. Anal. 52, 156 (2014).
  * Lucarini (2009) V. Lucarini, Journal of Statistical Physics 134, 381 (2009).
  * Morales  _et al._ (1999) C. A. Morales, M. J. Pacífico, and E. R. Pujals, Proc. Amer. Math. Soc. 127, 3393 (1999).
  * Reick (2002) C. H. Reick, Phys. Rev. E 66, 570 (2002).
  * Farmer (1985) D. J. Farmer, Phys. Rev. Letters 55, 351 (1985).
  * Ershov (1993) S. V. Ershov, Physics Letters A 177, 180 (1993).
  * Blonigan and Wang (2018) P. J. Blonigan and Q. Wang, Journal of Computational Physics 354, 447 (2018).
  * Albers  _et al._ (2006) D. J. Albers, J. C. Sprott, and J. P. Crutchfield, Phys. Rev. E 74 (2006).
  * Cvitanović  _et al._ (2010) P. Cvitanović, R. L. Davidchack, and E. Siminos, SIAM Journal on Applied Dynamical Systems 9, 1 (2010).
  * Bojanczyk  _et al._ (1992) A. W. Bojanczyk, G. H. Golub, and P. Van Dooren, in  _Proceedings of SPIE - The International Society for Optical Engineering_ (Cornell University, Ithaca, United States, 1992) pp. 31–42.
  * Hench and Laub (1994) J. J. Hench and A. J. Laub, IEEE Trans. Automat. Contr. 39, 1197 (1994).
  * Golub and Van Loan (2012) G. H. Golub and C. F. Van Loan,  _Matrix Computations_ (JHU Press, 2012).
  * Kressner (2006) D. Kressner,  _Numerical Methods for General and Structured Eigenvalue Problems_ (Springer Science & Business Media, 2006).
  * Saad (2011) Y. Saad,  _Numerical Methods for Large Eigenvalue Problems_ , Revised Edition (SIAM, 2011).
  * Goldhirsch  _et al._ (1987) I. Goldhirsch, P.-L. Sulem, and S. A. Orszag, Physica D 27, 311 (1987).


# MapVisual: A Visualization Tool for Memory Access Patterns

Pavlos Aimoniotis, Maria-Rafaela Gkeka and Nikolaos Bellas Department of
Electrical and Computer Engineering
University of Thessaly
{paimoniotis, margkeka<EMAIL_ADDRESS>Pavlos Aimoniotis, Maria Rafaela Gkeka
and Nikolaos Bellas∗ 1{paimoniotis, margkeka<EMAIL_ADDRESS>
###### Abstract

Memory bandwidth is strongly correlated to the complexity of the memory access
pattern of a running application. To improve memory performance of
applications with irregular and/or unpredictable memory patterns, we need
tools to analyze these patterns during application development. In this work,
we present a software tool for the analysis and visualization of memory access
patterns. We perform memory tracing and profiling, we do data processing and
filtering, and we use visualization algorithms to produce three dimensional
graphs that describe the patterns both in space and in time. Finally, we
evaluate our toolflow on a variety of applications.

##  I INTRODUCTION

Memory has traditionally been a major performance bottleneck in computing
systems. Complex data structures and dynamic memory allocations may result in
complex access patterns scattered all over the memory space. Additionally, the
memory has a predefined size that does not allow an unlimited number of
programs to run concurrently. For many years, the "Memory Wall" has been a
major problem referring to the growing speed gap between CPU and DRAM. [1],
[2], [3] have mentioned that one way to improve memory performance is by
optimizing memory access patterns (MAP). The difference between an L1 cache-
hit, and a full miss resulting in main-memory access, is about 50ns to 100ns.
If algorithms make random walks, they are less likely to benefit from the
hardware support that hides this latency.

Using code analysis and visualization, our contribution provides new insights
to assist code developers to better understand memory access patterns. There
have been several tools for code and memory profiling, such as Intel Advisor
and gprof. Although these tools provide excellent profiling, they lack
visualization. The MapVisual implementation flow is described briefly in the
following paragraph.

First, (1) we use a memory trace file and we analyze it through the Gleipnir
tool. Second, (2) we parse the Gleipnir output file and (3) we go through a
data filtering process to produce analytics. Finally, (4) we perform
visualization techniques on the output of step (3).

##  II BACKGROUND

###  II-A Valgrind Framework

Valgrind is an open-source dynamic binary instrumentation (DBI) framework used
for memory debugging, memory leak detection and profiling. Valgrind is a
Virtual Machine that uses Just-In-Time compilation techniques to translate the
assembly code into a processor-agnostic Intermediate Representation (IR). A
Valgrind tool (such as Lackey) can process the IR before Valgrind translates
the IR back into machine code for host-based execution.

###  II-B Gleipnir

Gleipnir is a memory analysis tool that maps an application’s source code
variables to generated data traces [4]. It is an extension of Valgrind’s
Lackey. The trace file generated by Lackey is purely a data and instruction
trace which includes only load and store instructions and their corresponding
data address. Gleipnir extends this information with the thread id,
originating function, scope, data structure and the corresponding element, or
a single scalar variable name. Gleipnir produces the following information
Operation, Address, Memory_Size, Thread_Id, Scope, Function, Variable_Info.

Gleipnir traces and produces an output line for every load, store and modify
in memory. This tracing can lead to a significant run time overhead. MapVisual
optimizes the Gleipnir tool specifically for our purposes by filtering the
amount of information it produces.

Ginwidth=0.3

(a) Blocked Matrix Multiplication (BMM) Complete Memory Access Pattern
Visualization

(b) Blocked MM C = A*B, Array B 2D Visualization

(c) 3D Array Accessed Sequentially Visualization

Figure 1: Memory Access Patterns Visualization for three Applications

##  III MapVisual Analysis

###  III-A Memory Tracing

Memory tracing is crucial for the performance of the tool. We run Gleipnir to
generate the information on the allocated memory, the variable name and
address starting point. We used the following two flags:

⬇

–read-var-info=yes

–read-debug=yes

However, this creates major bottlenecks regarding time and space performance.
We reduce the generated trace file $25\%$ on average, with bigger applications
tending to be around $75\%$. We observed a peak reduction of $83\%$ on one of
our benchmarks, achieving significant timing performance speedup.

###  III-B Data Structure

We parse the output file of Gleipnir in order to minimize the amount of
information produced [5]. We create record instances for every valid address
from the trace file and a Look-Up Table (LUT) for every unique record. For
every duplicate of an address, we create a list with the (1) Operation Type,
representing stores, loads and modifies, and the (2) LUT ID, representing the
position of that access in the LUT to avoid performance overhead of storing
the same data for every address.

###  III-C Data Filtering

The goal is to provide analytics prior to visualization. We determine for
individual memory addresses Memory Usage, Appearances and Variable Name.

Memory usage refers to statistics collected about operations per address.
Appearances refers to how many times each address has been accessed during the
execution of the program. Variable Name refers to the name of an array. The
process of extracting the name is mandatory for heap allocated elements, as
Gleipnir does not provide such information straightforward.

We memory traced and stored data sequentially, as a consequence tool
structures represent real-time accessibility. A timeline of the execution is
generated using this particular information.

###  III-D Visualization

We produce three types of visualization.

  1. 1.

Complete memory access pattern, refers to the whole execution of the program
and is the standard visualization of the tool. This 3D visualization
demonstrates for each variable the addresses and their usage.

  2. 2.

Two-dimension array visualization, refers to 2D arrays. It reconstructs an
array in a three-dimensional world. The third dimension represents information
about memory usage.

  3. 3.

Three-dimension array visualization, refers to 3D arrays. It reconstructs an
array in a three-dimensional world, where the axes show the position of each
element.

To visualize, we process the data we created of ordered information in III-C
and allocated information in III-B. Visualizations have different coloring, a
heat map that refers to the execution timeline of the program. Also, multiple
features are supported, such as zoom in, zoom out, and mouse-over information.

##  IV EXPERIMENTAL RESULTS

Figure 1 depicts the visualizations of three example benchmarks. Each one
corresponds to the types mentioned in III-D. The color heat map shows how far
in the past a specific address was accessed: a darker color indicates an older
access, whereas a brighter color indicates a more recent memory access.

Figure 1a shows the Complete MAP of a Blocked Matrix Multiplication (BMM)
$C=A*B$, axes $z$ and $y$ represent addresses and variables respectively. The
$x-axis$ (depth) represents the number of accesses per address and the color
represents the occurrence time of an access during the execution of the
program. Using this information we can observe the blocks accessed during the
multiplication. Figure 1b shows a reconstructed 2D array of BMM and how this
array has been accessed. In this visualization, axes represent rows $(z-axis)$
and columns $(x-axis)$ of an array, and number of appearances per array
element $(y-axis)$. In particular, this visualization represents array B of
$C=A*B$ BMM, and shows the memory access pattern in this particular array.
From the different coloring and number of accesses per element we can
determine the different pattern and how code utilized this array. Element
(0,0) that represents the start of the array is located at the bottom left,
due to Plotly library implementation. Finally, Figure 1c shows a
reconstruction of a 3D array that has been accessed sequentially. Axes
$(z,x,y)$ represent the coordinates of 3D array elements. This particular
array has been accessed first in depth $y-axis$, then towards columns $x-axis$
and finally towards rows $z-axis$. This is inferred from the slow color change
along the depth axis which means those elements are accessed near in time, the
slightly more frequent color change in columns, and finally the abrupt color
change in rows, where we conclude that there has been a long interval between
the two accesses.

We also evaluated our work in various other applications, such as Raycast
Algorithm and Mandelbrot Set.

##  V CONCLUSIONS

In this work, we presented a software for the visualization of Memory Access
Patterns. We demonstrated the memory tracing process, the implementation of
the tool and ideas on visualizing memory access patterns. Although further
optimizations can be applied, we provide a strong tool that can be crucial on
MAP analysis.

## References

  * [1] P. Grun, N. Dutt, N. Dutt, N. Dutt, and A. Nicolau, “Apex: access pattern based memory architecture exploration,” in _Proceedings of the 14th international symposium on Systems synthesis_. ACM, 2001, pp. 25–32.
  * [2] P. A. Boncz, S. Manegold, M. L. Kersten _et al._ , “Database architecture optimized for the new bottleneck: Memory access,” in _VLDB_ , vol. 99, 1999, pp. 54–65.
  * [3] B. Jang, D. Schaa, P. Mistry, and D. Kaeli, “Exploiting memory access patterns to improve memory performance in data-parallel architectures,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 22, no. 1, pp. 105–118, 2010.
  * [4] T. Janjusic and K. Kavi, “Gleipnir: A memory profiling and tracing tool,” _ACM SIGARCH Computer Architecture News_ , vol. 41, no. 4, pp. 8–12, 2013\.
  * [5] C. Ntogkas, “Design and implementation of a tool for memory access pattern visualization,” _University of Thessaly, Diploma Thesis_ , 2017. [Online]. Available: https://www.e-ce.uth.gr/wp-content/uploads/formidable/59/Ntogkas_Christos.pdf