k})]+Pr(\overline{E})\mathbb{E}[(\operatorname{Rew}(\pi^{\operatorname{OPT}}_{T,\infty})-\operatorname{Rew}(\pi^{\operatorname{OPT}}_{T,\leq
k}))|\overline{E}]\\\
&\leq\mathbb{E}[\operatorname{Rew}(\pi^{\operatorname{OPT}}_{T,k})]+Pr(\overline{E})\\\
&\leq\mathbb{E}[\operatorname{Rew}(\pi^{\operatorname{OPT}}_{T,k})]+\alpha^{k}\end{split}$

The last inequality follows from Lemma 3.

∎

###### Proof of Theorem 4.

For $\alpha\leq 0.29$, we have

$0\leq\operatorname{Rew}(\pi^{\operatorname{OPT}})-\operatorname{Rew}(\pi^{\operatorname{OPT}}_{T,\infty})\leq\alpha^{2}\cdot\frac{2-2\alpha+\alpha^{2}}{1-\alpha+\alpha^{2}}\cdot[\alpha\cdot\frac{2-\alpha}{1-\alpha}]^{T-2}$

and

$0\leq\mathbb{E}[\operatorname{Rew}(\pi^{\operatorname{OPT}}_{T,\infty})]-\mathbb{E}[\operatorname{Rew}(\pi^{\operatorname{OPT}}_{T,k})]\leq\alpha^{k}$

The theorem follows by adding the two inequalities. ∎

###  E.3. Locating the Root of Equation 2

In this section, we will approximate the root of the equation

$\begin{split}\nu&=\alpha(1-\alpha)+(1-\alpha)\nu^{2}\sum_{i=0}^{k-1}\nu^{i}\\\
&=\alpha(1-\alpha)+(1-\alpha)\nu^{2}\frac{1-\nu^{k}}{1-\nu}\end{split}$

The left hand side is the equation of a straight line. The right hand side of
the first equality is clearly convex in $\nu$. A straight line can intersect a
convex curve at most twice. Thus, the above equation can have at most two
roots.

At $\nu=0$, the LHS is zero and thus, is smaller than the RHS, equal to
$\alpha(1-\alpha)$. At $\nu=\infty$, the LHS is smaller than the RHS since the
RHS is a degree $k+1$ polynomial with a positive leading coefficient while the
LHS is just a line. At $\nu=\alpha$, the LHS equals $\alpha$. The RHS equals

$\alpha(1-\alpha)+(1-\alpha)\alpha^{2}\frac{1-\alpha^{k}}{1-\alpha}\leq\alpha(1-\alpha)+\alpha^{2}=\alpha$

and thus is smaller than the LHS. Thus, one of the roots of the equation is in
$[0,\alpha]$ and the other is in $[\alpha,\infty]$. Remember
$\nu=\alpha\,e\leq\alpha\times 1$, since $e$ is the probability of $E$
happening. Thus, we are interested in locating the smaller root of the above
equation.

###### Lemma 0.

For $\alpha\leq 0.29$, the smaller root of

(3) $\nu=\alpha(1-\alpha)+(1-\alpha)\nu^{2}\frac{1-\nu^{k}}{1-\nu}$

is at least $\alpha-\alpha^{k+1}$.

###### Proof.

Consider the equation

(4) $\nu=\alpha(1-\alpha)+(1-\alpha)\nu^{2}\frac{1-\alpha^{k}}{1-\nu}$

instead. The right hand side is still convex in $[0,1]$. In the region
$[0,\alpha]$, the RHS of the above equation is strictly smaller than
$\alpha(1-\alpha)+(1-\alpha)\nu^{2}\frac{1-\nu^{k}}{1-\nu}$. Thus, at the root
$\nu_{1}$ of Equation 3

$\nu_{1}>\alpha(1-\alpha)+(1-\alpha)\nu_{1}^{2}\frac{1-\alpha_{1}^{k}}{1-\nu_{1}}$

The left hand side of the new equation has already overtaken the right hand
side at $\nu=\nu_{1}$. The root of Equation 4 is smaller than $\nu_{1}$. We
will find this root and use it to lower bound $\nu_{1}$.

Next, at $\nu=\alpha-\alpha^{k+1}$, we will show that the RHS of the new
equation is larger than the LHS. Since the LHS overtakes the RHS at
$\nu=\alpha$, the root will have to lie in $[\alpha-\alpha^{k+1},\alpha]$. For
$\alpha\leq 0.29$, the following the chain of inequalities hold.

$\begin{split}0&\geq-2\alpha^{2}+4\alpha-1\\\
(1-\alpha)&\geq\alpha^{2}+3(1-\alpha)\alpha\\\
(1-\alpha)\alpha^{k-1}&\geq\alpha^{k+1}+3(1-\alpha)\alpha^{k}\\\
(1-\alpha)\alpha^{k-1}&\geq\alpha^{k+1}+(1-\alpha)\alpha^{k}[1+(1-\alpha^{k})+(1-\alpha^{k})^{2}]\\\
(1-\alpha+\alpha^{k+1})\times\alpha^{k-1}&\geq(1-\alpha)(1-[1-\alpha^{k}]^{3})+\alpha^{k+1}\\\
\alpha^{k+1}&\geq\alpha^{2}\times[1-(1-\alpha)\frac{(1-\alpha^{k})^{3}}{1-\alpha+\alpha^{k+1}}]\end{split}$

Showing

$\alpha-\alpha^{k+1}\leq\alpha-\alpha^{2}\times[1-(1-\alpha)\frac{(1-\alpha^{k})^{3}}{1-\alpha+\alpha^{k+1}}]=\alpha(1-\alpha)+(1-\alpha)[\alpha-\alpha^{k+1}]^{2}\times\frac{1-\alpha^{k}}{1-\alpha+\alpha^{k+1}}$

is now straightforward, completing the proof. ∎

##  Appendix F Errors in Estimation through Finite Sampling

Given a distribution $\operatorname{D}_{0}$, we want to approximate the
sequence of distributions $\big{(}\operatorname{D}_{t}\big{)}_{0\leq t\leq T}$
such that
$\operatorname{D}_{t+1}=\text{AddLayer}(\alpha,\beta,\lambda,k,\pi^{\operatorname{OPT}}_{T,k},\operatorname{D}_{t})$.

###  F.1. Bounding Expected Rewards through Chernoff Bound and the McDiarmid’s
Inequality

Note that we need an estimate of the CDF of each $\operatorname{D}_{t}$ and
just getting an estimate of the expected value of $\operatorname{D}_{t}$ is
not sufficient. While the Chernoff bound only gives an estimate of the
expected reward of $\operatorname{D}_{t+1}$ given the distribution
$\operatorname{D}_{t}$, using the McDiarmid’s inequality to directly get tail
bounds on the expected value of $\operatorname{D}_{T}$ requires an exponential
number of samples in $T$ as discussed below.

As a thought experiment, we use $n$ samples to get an approximation
$\hat{\operatorname{D}}_{1}$ of the distribution $\operatorname{D}_{1}$, and
allow
$\text{AddLayer}(\alpha,\beta,\lambda,k,\pi^{\operatorname{OPT}}_{T,k},\cdot)$
to construct infinitely many samples to approximate distributions from there
on to get
$\hat{\operatorname{D}}_{2},\hat{\operatorname{D}}_{3},\dots,\hat{\operatorname{D}}_{T-1},\hat{\operatorname{D}}_{T}$.
We sketch a ‘fake’ proof using the McDiarmid’s inequality to get an idea on
the number of samples that would be recommended by McDiarmid’s inequality if
we want to bound the error in estimating the expected reward by (say) $0.001$
with probability at least $99\%$.

###### Theorem 1 (McDiarmid’s inequality, McDiarmid et al.,, 1989).

Let
$f:\mathcal{X}_{1}\times\mathcal{X}_{2}\times\dot{\times}\mathcal{X}_{n}\xrightarrow{}\mathbb{R}$
satisfy the $(c_{1},\dots,c_{n})$-bounded difference property for
$\delta_{1},\delta_{2},\dots,\delta_{n}\geq 0$, i.e.,

$\sup_{x_{i}^{\prime}\in\mathcal{X}_{i}}|f(x_{1},x_{2},\dots,x_{i},\dots,x_{n})-f(x_{1},x_{2},\dots,x_{i}^{\prime},\dots,x_{n})|\leq\delta_{i}$

for all
$(x_{1},\dots,x_{n})\in\mathcal{X}_{1}\times\mathcal{X}_{2}\times\dot{\times}\mathcal{X}_{n}$
and for all $1\leq i\leq n$. Then, for independent random variables
$X_{1}\in\mathcal{X}_{1},X_{2}\in\mathcal{X}_{2},\dots,X_{n}\in\mathcal{X}_{n}$,
and $\epsilon\geq 0$

$Pr(|f(X_{1},\dots,X_{n})-\mathbb{E}[f(X_{1},\dots,X_{n})]|\geq\epsilon)\leq
2e^{-\tfrac{2\epsilon^{2}}{\sum_{i=1}^{n}\delta_{i}^{2}}}$

The entire simulation can be thought of as a function of the samples
$\big{(}s_{\ell}\big{)}_{1\leq\ell\leq n}$ drawn to construct
$\hat{\operatorname{D}}_{1}$. To use the McDiarmid’s inequality We want to
find a constant $\delta$ such that changing a single sample $s_{\ell}$ would
not change the outcome of the simulation by more than $\delta$.

Suppose we ‘taint’ a sample $s_{\ell}$ in the approximation
$\hat{\operatorname{D}}_{1}$ (this amounts to tainting $\tfrac{1}{n}$ fraction
of the samples constructed for $\hat{\operatorname{D}}_{1}$). For each of the
infinitely many samples we construct to get $\hat{\operatorname{D}}_{2}$ from
$\hat{\operatorname{D}}_{1}$, we draw a sample for each of the adversary’s
wallets with the $k$ smallest scores. Thus, $\tfrac{k}{n}$ fractions of the
samples constructed for $\hat{\operatorname{D}}_{2}$ would have used
$s_{\ell}$ and are tainted. This is ignoring the fact that we take an
expectation with respect $r_{0}$, the reward from remaining silent; the number
of tainted samples only gets worse if this gets accounted. Inductively,
$\frac{k^{T-1}}{n}$ samples of $\hat{\operatorname{D}}_{T}$ are tainted. All
tainted samples in round $T$ can potentially increase from $-T\,\lambda$ to
$T\,(1-\lambda)$. Thus, the expected value of the distribution can increase by
$\delta=\tfrac{k^{T-1}}{n}\,T$.

For $\delta$ described above, $T=5$ and $k=5$, and for an error in the
estimated reward of at most $0.001$ with probability $99\%$, we would require
$n=0.75\times 10^{11}$. The choices of $T,k$ and the confidence probability
are very mild and we need a huge number of samples just to account for the
error from finitely sampling $\hat{\operatorname{D}}_{1}$. The number of
samples would only grow if
$\hat{\operatorname{D}}_{2},\dots,\hat{\operatorname{D}}_{T}$ are also
constructed through finite sampling. Of course, the analysis is extremely
loose and can be improved. The number of samples recommended by any such
concentration bounds turn out to be extremely large.

###  F.2. Bias

We revisit the computation of the adversary’s reward when its $k$ smallest
scores are $\vec{c}_{-0}$ and the corresponding rewards are $\vec{r}_{-0}$.
The reward from broadcasting the credential $i$ (or remaining silent) equals
$e^{-c_{i}\,(1-\beta)\,(1-\alpha)}\,(r_{i}+\mathbbm{1}(i\neq 0))$ (assuming a
stopping time is not reached). The adversary chooses the action that maximizes
its reward.

Consider a toy version of the adversarial game over just two rounds. In the
first round, a coin with heads probability $\frac{1}{k}$ is tossed $n$ times.
Let the empirically observed probability of heads be $p_{1}$. For each trial
in the second round, toss $k$ coins each with heads probability $p_{1}$. The
outcome is $1$ even if one of $k$ tosses turns out heads (we take the maximum
amongst $k$ Bernoulli trials, similar to the adversarial game). Repeat the
trial $n$ times to observe an empirical probability $p_{2}$.

We verify that $p_{2}$ does not remain constant for $k=2$ over different
choices of $n$. If $n=\infty$, $p_{1}=\tfrac{1}{k}=\tfrac{1}{2}$ and
$p_{2}=1-(1-\tfrac{1}{k})^{2}=\tfrac{3}{4}$. However, for $n=10$ samples,
$p_{1}=\tfrac{m}{n}$ with probability $\binom{n}{m}\cdot 2^{-n}$ and
$\mathbb{E}[p_{2}]=\sum_{m=0}^{n}\big{[}1-\big{(}1-\tfrac{m}{n}\big{)}^{2}\big{]}\times\binom{n}{m}\cdot
2^{-n}\approx 0.784\neq 0.75$

A similar bias creeps into the adversarial game. The effect of the bias on the
estimated rewards worsens with the number of rounds $T$ and becomes better
with the number of samples $n$.

###  F.3. Proof of Theorem 5

We first show the following concentration inequality between the CDFs of a
distribution $D$ and a distribution constructed empirically by sampling from
$D$.

###### Theorem 2.

Let $\operatorname{D}$ be some distribution supported on
$[-\lambda,t\,(1-\lambda)]$ and let $\underline{\operatorname{D}}$ and
$\overline{\operatorname{D}}$ be the result of
$\text{Deflate}(n,\gamma,\hat{\operatorname{D}},t)$ and
$\text{Inflate}(n,\gamma,\omega,\hat{\operatorname{D}},t)$ respectively, where
$\hat{D}$ is constructed by sampling $n$ times from $\operatorname{D}$.

  1. (1)

$\underline{\operatorname{D}}$ is dominated by $\operatorname{D}$ with
probability at least $1-\gamma$.

  2. (2)

$\overline{\operatorname{D}}$ dominates $\operatorname{D}$ with probability at
least $1-\Big{(}\gamma+\tfrac{e^{-\omega
n}}{\omega}\,\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}\Big{)}$

For a distribution $\operatorname{D}$ and a quantile $q$, let the inverse CDF
$\operatorname{D}^{-1}(q)$ denote the value corresponding to the quantile
$q$121212Technically, the value corresponding to a quantile $q$ is given by
$D^{-1}(1-q)$. For cleanliness, we pretend that the CDF $D(s)$ describes the
probability $Pr_{S\sim D}(S\geq s)$ and thus, $D^{-1}(q)$ is the value
associated to the quantile $q$.

###### Theorem 3 (DKW inequality; Dvoretzky et al.,, 1956; Massart,, 1990).

Let $\operatorname{D}$ be some distribution and let $\hat{\operatorname{D}}$
be constructed by drawing $n$ samples from $\operatorname{D}$. Then,

$Pr\big{(}\hat{\operatorname{D}}^{-1}(q+\sqrt{\tfrac{\ln\gamma^{-1}}{2n}})\leq\operatorname{D}^{-1}(q)\leq\hat{\operatorname{D}}^{-1}(q-\sqrt{\tfrac{\ln\gamma^{-1}}{2n}})\text{
for all }q\in[0,1]\big{)}\geq 1-\gamma$

###### Lemma 0.

Let $\underline{\operatorname{D}}=\text{Deflate}(n,\gamma,\operatorname{D},t)$
for some distribution $\operatorname{D}$ with infimum $\psi$.
$\underline{\operatorname{D}}$ is dominated by $\operatorname{D}$ with
probability at least $1-\gamma$.

###### Proof.

Consider the distribution $\hat{\operatorname{D}}$ constructed from drawing
$n$ samples from $\operatorname{D}$. By the DKW inequality,
$\hat{\operatorname{D}}^{-1}(q+\sqrt{\tfrac{\ln\gamma}{2n}})\leq\operatorname{D}^{-1}(q)$.
By deleting the $\sqrt{\tfrac{\ln\gamma}{2n}}$ smallest (strongest) quantiles
of $\hat{\operatorname{D}}$ and inserting the infimum $-\lambda$ of
$\operatorname{D}$ to construct $\underline{\operatorname{D}}$, we decrease
the quantile of each sample by $\sqrt{\tfrac{\ln\gamma}{2n}}$. Thus,

$\underline{\operatorname{D}}^{-1}(q)=\hat{\operatorname{D}}^{-1}(q++\sqrt{\tfrac{\ln\gamma}{2n}})\leq\operatorname{D}^{-1}(q)$

for all $q\in[0,1-\sqrt{\tfrac{\ln\gamma}{2n}}]$. For
$q\in[1-\sqrt{\tfrac{\ln\gamma}{2n}},1]$,

$\underline{\operatorname{D}}^{-1}(q)=\psi\leq\operatorname{D}^{-1}(q)$

∎

###### Lemma 0.

Let
$\overline{\operatorname{D}}=\text{Inflate}(n,\gamma,\omega,\operatorname{D},t)$
for some distribution $\operatorname{D}$ with supremum $\Psi$.
$\overline{\operatorname{D}}$ dominates $\operatorname{D}$ with probability at
least $1-\Big{(}\gamma+\tfrac{e^{-\omega
n}}{\omega}\,\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}\Big{)}$.

###### Proof.

Let $\hat{\operatorname{D}}$ be the distribution constructed by sampling $n$
times from $\operatorname{D}$. We show that, with very high probability, a
sample with quantile $q\in(m\,\omega,(m+1)\,\omega]$ was drawn while
constructing $\hat{\operatorname{D}}$ for all $0\leq
m<\frac{1}{\omega}\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$. Inserting $\omega
n$ of each such sample with extremely small quantiles would ensure
$\overline{\operatorname{D}}$ dominates $\operatorname{D}$.

Consider a single draw from $\operatorname{D}$. The probability that the
outcome of the draw does not have a quantile $q\in(m\,\omega,(m+1)\,\omega]$
equals $1-\omega$. The probability that none of the $n$ samples have a
quantile in the range $(m\,\omega,(m+1)\,\omega]$ equals $(1-\omega)^{n}\leq
e^{-\omega\,n}$. By a union bound, the probability that there does not exist
any sample with quantile $q\in(m\,\omega,(m+1)\,\omega]$ for some $0\leq
m<\frac{1}{\omega}\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$ is at most
$e^{-\omega n}\times\frac{1}{\omega}\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$.

Next, we show that
$\overline{\operatorname{D}}^{-1}(q)\geq\operatorname{D}^{-1}(q)$ for
$q\in[0,\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}]$ with probability at least
$1-\tfrac{e^{-\omega n}}{\omega}\,\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$. As
argued above, with probability at least $1-\tfrac{e^{-\omega
n}}{\omega}\,\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$, the $m\textsuperscript{th}$
largest sample $s_{m}$ of $\hat{\operatorname{D}}$ has a quantile $q\leq
m\omega$ for $1\leq m<\frac{1}{\omega}\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$.
We delete the $n\,\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$ smallest samples and
append $\omega n$ copies each of $s_{0}:=\Psi,s_{1},\dots,s_{M}$ for
$M=\Big{(}\frac{1}{\omega}\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}-1\Big{)}$.
For a quantile $q\in(m\omega,(m+1)\,\omega]$ ($0\leq
m\leq\frac{1}{\omega}\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$),
$\overline{\operatorname{D}}^{-1}(q)=s_{m}\geq\operatorname{D}^{-1}(q)$.

For $q\in(\sqrt{\tfrac{\ln\gamma^{-1}}{2n}},1]$, with probability at least
$1-\gamma$,
$\hat{\operatorname{D}}^{-1}\big{(}q-\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}\big{)}\geq\operatorname{D}^{-1}(q)$
(by the DKW inequality). We appended $n\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$
samples to $\hat{\operatorname{D}}$, all of them with quantiles at most
$\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$, to get $\overline{\operatorname{D}}$.
Thus,
$\overline{\operatorname{D}}^{-1}(q)=\hat{\operatorname{D}}^{-1}(q-\sqrt{\tfrac{\ln\gamma^{-1}}{2n}})\geq\operatorname{D}^{-1}(q)$.

Combining the two, we get that the probability that
$\overline{\operatorname{D}}$ dominates $\operatorname{D}$ is at least
$1-\Big{(}\gamma+\tfrac{e^{-\omega
n}}{\omega}\,\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}\Big{)}$. ∎

Lemma 4 and Lemma 5 prove Theorem 2.

Thus, to prove Theorem 5, it is sufficient to show that the reward
distributions are always within the range $[-\lambda,t\,(1-\lambda)]$.

###### Proof of Theorem 5.

Clearly, the maximum reward the adversary can earn in each round is
$(1-\lambda)$ and the supremum of
$\operatorname{D}^{\operatorname{OPT}}_{t,k}$ for each round $t$ is at most
$t\,(1-\lambda)$. Combining Theorem 2 with a straightforward union bound
between the $T$ inflation procedures performed across the $T$ rounds, we get
the first part.

For the second part, a similar proof to the first will follow provided we show
the infimum of $\hat{\operatorname{D}}^{\operatorname{OPT}}_{t,k}$ is
$-\lambda$ for each $1\leq t\leq T$. We show this by sketching a strategy that
guarantees a reward at least $-\lambda$ for the adversary.
$\pi^{\operatorname{OPT}}_{T,k}$ being the optimal strategy will only provide
the adversary with a better reward. Consider the strategy where the adversary
keeps broadcasting the credential of a random adversarial potential winner
until a stopping time is reached. The adversary wins $(1-\lambda)$ each round
before the first stopping time and $-\lambda$ for the stopping time, since the
adversary plays the round without getting elected. Since $\lambda<1$, the
adversary’s reward is at least $-\lambda$. Note that the adversary guarantees
itself a reward greater than $-\lambda$ even conditioned on the credentials of
the current and future rounds. Thus, the optimal strategy
$\pi^{\operatorname{OPT}}_{T,k}$ can adapt the above strategy if the adversary
faces an extremely unfavourable seed and achieve a reward at least $-\lambda$.

∎

##  Appendix G Computing Expectations

We get into the details next. The first step towards pre-computing is
chiselling the expression for the samples into a form that allows intensive
reuse of values.

$\begin{split}s_{\ell}=\mathbb{E}_{c_{0}\sim\operatorname{exp}(\beta\,(1-\alpha)),r_{0}\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\big{[}\max_{0\leq
i\leq
i^{*}(\vec{c})}\\{e^{-c_{i}\cdot(1-\beta)\cdot(1-\alpha)}\,(r_{i}+\mathbbm{1}(i\neq
0))\\}\cdot\mathbbm{1}(i^{*}(\vec{c})\neq 0)\big{]}-\lambda\end{split}$

Remember that $i^{*}(\vec{c})$ is the number of adversarial potential winners,
$h(c_{0},r_{0})=e^{-c_{0}\cdot(1-\beta)(1-\alpha)}r_{0}$ is the reward the
adversary gets by staying silent and
$g(c_{0},\vec{c}_{-0},\vec{r}_{-0})=\max_{0\leq i\leq
i^{*}(\vec{c})}{e^{-c_{i}\cdot(1-\beta)(1-\alpha)}}(1+r_{i})$ is the maximum
reward the adversary gets by broadcasting one of its credentials. Thus,

$\begin{split}s_{\ell}&=\mathbb{E}_{c_{0}\leftarrow\operatorname{exp}(\beta(1-\alpha)),r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(c_{0},\vec{c}_{-0},\vec{r}_{-0})\\}\mathbbm{1}\\{i^{*}(\vec{c})\neq
0\\}\right]-\lambda\\\
&=\int_{0}^{\infty}\mathbb{E}_{r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(c_{0},\vec{c}_{-0},\vec{r}_{-0})\\}\mathbbm{1}\\{i^{*}(\vec{c})\neq
0\\}\right]\beta\,(1-\alpha)\,e^{-c_{0}\,\beta\,(1-\alpha)}\,dc_{0}-\lambda\\\
\end{split}$

As a first step, we break down the integral wrt $c_{0}$ as a sum of integrals
over smaller intervals where $i^{*}(\vec{c})$ remains constant, i.e.,

$\begin{split}s_{\ell}=\int_{0}^{c_{1}}&\mathbb{E}_{r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(c_{0},\vec{c}_{-0},\vec{r}_{-0})\\}\mathbbm{1}\\{i^{*}(\vec{c})\neq
0\\}\right]\beta\,(1-\alpha)\,e^{-c_{0}\,\beta\,(1-\alpha)}\,dc_{0}\\\
&+\sum_{i^{*}=1}^{k}\int_{c_{i^{*}}}^{c_{{i^{*}}+1}}\mathbb{E}_{r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(c_{0},\vec{c}_{-0},\vec{r}_{-0})\\}\mathbbm{1}\\{i^{*}(\vec{c})\neq
0\\}\right]\beta\,(1-\alpha)\,e^{-c_{0}\,\beta\,(1-\alpha)}\,dc_{0}-\lambda\end{split}$

Observe that when $c_{0}\in[c_{i^{*}},c_{i^{*}+1})$, $i^{*}(\vec{c})=i^{*}$.
In this range,

$\begin{split}g(c_{0},\vec{c}_{-0},\vec{r}_{-0})&=\max_{1\leq i\leq
i^{*}(\vec{c})}{e^{-c_{i}\cdot(1-\beta)(1-\alpha)}}(1+r_{i})\\\ &=\max_{1\leq
i\leq i^{*}}{e^{-c_{i}\cdot(1-\beta)(1-\alpha)}}(1+r_{i})\end{split}$

and is independent of the value of $c_{0}$. When $c_{0}\in[0,c_{1})$,
$i^{*}(\vec{c})=0$ and so is
$\mathbb{E}_{r_{0}\leftarrow
D}\left[\max\\{h(c_{0},r_{0}),g(c_{0},\vec{c}_{-0},\vec{r}_{-0})\\}\mathbbm{1}\\{i^{*}(\vec{c})\neq
0\\}\right]$. The entire integral wrt $c_{0}$ over $[0,c_{1})$ equals zero.
$s_{\ell}$ simplifies to

$s_{\ell}=\sum_{i^{*}=1}^{k}\int_{c_{i^{*}}}^{c_{{i^{*}}+1}}\mathbb{E}_{r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(c_{0},\vec{c}_{-0},\vec{r}_{-0})\\}\right]\beta\,(1-\alpha)\,e^{-c_{0}\,\beta\,(1-\alpha)}\,dc_{0}-\lambda$

Let
$f(i^{*},\vec{c}_{-0},\vec{r}_{-0})=\int_{c_{i^{*}}}^{c_{{i^{*}}+1}}\mathbb{E}_{r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(c_{0},\vec{c}_{-0},\vec{r}_{-0})\\}\right]\beta\,(1-\alpha)\,e^{-c_{0}\,\beta\,(1-\alpha)}\,dc_{0}$
and therefore,

$s_{\ell}=\sum_{i^{*}=1}^{k}f(i^{*},\vec{c}_{-0},\vec{r}_{-0})-\lambda$

We speed up the computation of $f(i^{*},\vec{c}_{-0},\vec{r}_{-0})$ next.
$f(i^{*},\vec{c}_{-0},\vec{r}_{-0})$ is computed by integrating $c_{0}$ over
$[c_{i^{*}},c_{i^{*}+1})$, where, as discussed in the previous paragraph, the
value of $g(c_{0},\vec{c}_{-0},\vec{r}_{-0})$ is independent of $c_{0}$. To
signify this independence, we will re-parameterize
$g(c_{0},\vec{c}_{-0},\vec{r}_{-0})$ into
$g(i^{*},\vec{c}_{-0},\vec{r}_{-0})$.

$g(i^{*},\vec{c}_{-0},\vec{r}_{-0})=\max_{1\leq i\leq
i^{*}}{e^{-c_{i}\cdot(1-\beta)(1-\alpha)}}(1+r_{i})$
$f(i^{*},\vec{c}_{-0},\vec{r}_{-0})=\int_{c_{i^{*}}}^{c_{{i^{*}}+1}}\mathbb{E}_{r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\\}\right]\beta\,(1-\alpha)\,e^{-c_{0}\,\beta\,(1-\alpha)}\,dc_{0}$

We set up pre-computations for
$\mathbb{E}_{r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\\}\right]$.
Note that we are computing
$\mathbb{E}_{r_{0}\leftarrow\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}\left[\max\\{h(c_{0},r_{0}),g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\\}\right]$
conditioned on $c_{0}$ and $i^{*}$.

$\begin{split}\mathbb{E}_{r_{0}\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}[\max\\{h(c_{0},r_{0}),g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\\}]&=g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\times
Pr(h(c_{0},r_{0})\leq g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))\\\ &\hskip
28.45274pt+\mathbb{E}_{r_{0}\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}[h(c_{0},r_{0})\times\mathbbm{1}(h(c_{0},r_{0})>g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))]\\\
&=g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\times
Pr(r_{0}\leq\tfrac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}\cdot(1-\beta)(1-\alpha)}})\\\
&\hskip
28.45274pt+\mathbb{E}_{r_{0}\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}[e^{-c_{0}\cdot(1-\beta)(1-\alpha)}r_{0}\times\mathbbm{1}(r_{0}>\tfrac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}\cdot(1-\beta)(1-\alpha)}})]\\\
&=g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\times
Pr(r_{0}\leq\tfrac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}\cdot(1-\beta)(1-\alpha)}})\\\
&\hskip
28.45274pt+e^{-c_{0}\cdot(1-\beta)(1-\alpha)}\times\mathbb{E}_{r_{0}\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}[r_{0}\times\mathbbm{1}(r_{0}>\tfrac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}\cdot(1-\beta)(1-\alpha)}})]\\\
&=e^{-c_{0}\cdot(1-\beta)(1-\alpha)}\times\Big{[}\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}\cdot(1-\beta)(1-\alpha)}}\times
Pr(r_{0}\leq\tfrac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}\cdot(1-\beta)(1-\alpha)}})\\\
&\hskip
28.45274pt+\mathbb{E}_{r_{0}\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}[r_{0}\times\mathbbm{1}(r_{0}>\tfrac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}\cdot(1-\beta)(1-\alpha)}})]\Big{]}\end{split}$

The term inside the big square in the final equation is of the form $\\{\theta
Pr(r_{0}\leq\theta)+\mathbb{E}_{r_{0}\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}[r_{0}\times\mathbbm{1}(r_{0}>\theta)]\\}$,
and can be pre-computed without knowing $\vec{c}_{-0}$ and $\vec{r}_{-0}$ (up
to a discretization error $\epsilon$). We produce a recipe to set-up pre-
computations for a general distribution $\operatorname{D}$.

$\text{Precompute}(\hat{D},\epsilon,\eta,t)$: (1) Construct the pdf $d$ of
$\hat{D}$ up to a discretization error $\epsilon$. (2) Construct the cdf
$\hat{D}$ recursively using the following recursion, once again, up to a
discretization error $\epsilon$.
$\hat{\operatorname{D}}(\theta-\epsilon)=Pr(r_{0}\leq\theta-\epsilon)=\hat{D}(\theta)-d(\theta)$
(3) Construct the right tail of the expectation
$E(\theta)=\mathbb{E}_{r_{0}\sim D}[r_{0}\times\mathbbm{1}(r_{0}>\theta)]$
recursively by $E(\theta-\epsilon)=E(\theta)+(\theta-\epsilon)\times
d(\theta-\epsilon)$ (4) Compute
$E_{\max}(\theta)=\theta\hat{\operatorname{D}}(\theta)+E(\theta)$ for all
$-t\lambda\leq\theta\leq t(1-\lambda)$ up to a discretization error $\epsilon$

The above compute takes time $O(\frac{t}{\epsilon})$ since the range of the
reward distribution has a width at $t$ in round $t$. Plugging this back into
the equation for $f(i^{*},\vec{c}_{-0},\vec{r}_{-0})$,

$\begin{split}f(i^{*},\vec{c}_{-0},\vec{r}_{-0})&=\int_{c_{i^{*}}}^{c_{i^{*}+1}}\mathbb{E}_{r_{0}\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}}[\max\\{h(c_{0},r_{0}),g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\\}]\beta\,(1-\alpha)\,e^{-c_{0}\,\beta\,(1-\alpha)}\,dc_{0}\\\
&=\int_{c_{i^{*}}}^{c_{i^{*}+1}}\beta(1-\alpha)e^{-c_{0}\beta(1-\alpha)}\times
e^{-c_{0}\,(1-\beta)\,(1-\alpha)}\times
E_{\max}(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}(1-\beta)(1-\alpha)}})\,dc_{0}\end{split}$

The final equality follows since $E_{\max}(\cdot)$ was defined such that
$\mathbb{E}_{r_{0}\sim
D}[\max\\{h(c_{0},r_{0}),g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\\}]=e^{-c_{0}\,(1-\beta)\,(1-\alpha)}\times
E_{\max}(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}(1-\beta)(1-\alpha)}})$.
We have all the ingredients to compute $f(i^{*},\vec{c}_{-0},\vec{r}_{-0})$
and speeding up the computation of the integral wrt $c_{0}$ is the final
frontier. At a high level, we compute

$G(c,g(i^{*},\vec{c}_{-0},\vec{r}_{-0})):=\int_{0}^{c}\beta(1-\alpha)e^{-c_{0}\beta(1-\alpha)}\times
e^{-c_{0}\,(1-\beta)\,(1-\alpha)}\times
E_{\max}(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}(1-\beta)(1-\alpha)}})\,dc_{0}$

and thus,

$f(i^{*},\vec{c}_{-0},\vec{r}_{-0})=G(c_{i^{*}+1},g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))-G(c_{i^{*}},g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))$

Case-1: $\beta=1$

This is the simplest case. The input to $E_{\max}(\cdot)$ becomes independent
of $c_{0}$.

$\begin{split}G(c,g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))&=\int_{0}^{c}(1-\alpha)e^{-c_{0}(1-\alpha)}\times
E_{\max}(g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))\,dc_{0}\\\
&=E_{\max}(g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))(1-e^{-c\,(1-\alpha)})\end{split}$

In this case, $f(i^{*},\vec{c}_{-0},\vec{r}_{-0})$ can be explicitly computed.

$f(i^{*},\vec{c}_{-0},\vec{r}_{-0})=E_{\max}(g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))\big{[}e^{-c_{i^{*}}\,(1-\alpha)}-e^{-c_{i^{*}+1}\,(1-\alpha)}\big{]}$

All of our pre-compute gets done in $O(\frac{t}{\epsilon})$ time.

Summarizing,

$\begin{split}s_{\ell}&=\sum_{i^{*}=1}^{k}f(i^{*},\vec{c}_{-0},\vec{r}_{-0})-\lambda\\\
&=\sum_{i^{*}=1}^{k}E_{\max}(g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))\big{[}e^{-c_{i^{*}}\,(1-\alpha)}-e^{-c_{i^{*}+1}\,(1-\alpha)}\big{]}-\lambda\end{split}$

Remember that $E_{\max}(\theta)=\theta D(\theta)+E(\theta)$ has already been
pre-computed.

Case-2: $\beta=0$

This is another straightforward case. When $\beta=0$, the adversary cannot
abstain from broadcasting. Thus, the adversarial reward equals
$g(k,\vec{c}_{-0},\vec{r}_{-0})$.

$s_{\ell}=g(k,\vec{c}_{-0},\vec{r}_{-0})-\lambda$

In particular, we circumvent computing $f(i^{*},\vec{c}_{-0},\vec{r}_{-0})$
altogether for $\beta=0$.

Case-3: $\beta\notin\\{1,0\\}$.

The most obvious way to compute $G(c,g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))$
would be to discretize the integral wrt to $c_{0}$ into steps of $\eta$.
However, the domain of $c_{0}$ is $[0,\infty)$.
$E_{\max}(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}(1-\beta)(1-\alpha)}})$
becomes constant for a sufficiently large value $C_{0}$ of $c_{0}$, but, it is
quite wasteful to compute the integral in discrete steps of $\eta$ from $0$ to
a very large $C_{0}$. Instead, we change variables and integrate wrt
$e^{-c_{0}\,\beta\,(1-\alpha)}$ and discretize this integral by $\eta$.

The following expressions just substantiate the above intuition.

$\begin{split}G(c,g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))&=\int_{0}^{c}\beta(1-\alpha)e^{-c_{0}\beta(1-\alpha)}\times
e^{-c_{0}\,(1-\beta)\,(1-\alpha)}\times
E_{\max}(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}(1-\beta)(1-\alpha)}})\,dc_{0}\\\
&=\int_{e^{-c\,\beta\,(1-\alpha)}}^{1}e^{-c_{0}\,(1-\beta)\,(1-\alpha)}E_{\max}(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{e^{-c_{0}(1-\beta)(1-\alpha)}})\,de^{-c_{0}\,\beta\,(1-\alpha)}\\\
&=\int_{e^{-c\,\beta\,(1-\alpha)}}^{1}\zeta^{\frac{1-\beta}{\beta}}E_{\max}\left(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{\zeta^{\frac{1-\beta}{\beta}}}\right)\,d\zeta\end{split}$

For every possible value
$g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\in[-t\,\lambda,t\,(1-\lambda)]$ and a
discretization parameter $\eta$, the above integral can be discretized as
follows.

$\begin{split}G(c,g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))=\eta\times\sum_{\zeta=0,\text{
in steps of
}\eta}^{c}\zeta^{\frac{1-\beta}{\beta}}E_{\max}\left(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{\zeta^{\frac{1-\beta}{\beta}}}\right)\end{split}$

Remember to take the upper and lower Riemann sums while computing
$\overline{\operatorname{D}}_{t}$ and $\underline{\operatorname{D}}_{t}$
respectively.

Once again, we compute the sample $s_{\ell}$ by

$\begin{split}s_{\ell}&=\sum_{i^{*}=1}^{k}G(c_{i^{*}+1},g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))-G(c_{i^{*}},g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))-\lambda\\\
&=\sum_{i^{*}=1}^{k}f(i^{*},\vec{c}_{-0},\vec{r}_{-0})-\lambda\end{split}$

The integrals with respect to $c_{0}$ take $O(\frac{1}{\eta})$ to pre-compute
and we do this once for each of the $\big{(}\frac{t}{\epsilon}\big{)}$
possible values of $g(i^{*},\vec{c}_{-0},\vec{r}_{-0})$, amounting to a total
run-time $O(\tfrac{t}{\epsilon\eta})$. Lemma 6 follows.

By Lemma 6, the pre-computations in one execution of
$\text{FiniteSampleAddLayer}(\alpha,\beta,\lambda,k,t,n,\pi^{\operatorname{OPT}}_{T,k},\gamma,\omega,\hat{\operatorname{D}}^{\operatorname{OPT}}_{t,k})$
takes $O(\tfrac{t}{\epsilon\,\eta})$ time. Across the $T$ rounds, a total
$O(\tfrac{T^{2}}{\epsilon\,\eta})$ time is spent in pre-computations.
Computing one sample $s_{\ell}$ involves summing up $k$ terms,

$s_{\ell}=\sum_{i^{*}=1}^{k}f(i^{*},\vec{c}_{-0},\vec{r}_{-0})-\lambda$

and takes $O(k)$ time to compute. A total of $O(kn)$ time is spent on sampling
in each execution of
$\text{FiniteSampleAddLayer}(\alpha,\beta,\lambda,k,t,n,\pi^{\operatorname{OPT}}_{T,k},\gamma,\omega,\cdot)$
and a total of $O(Tkn)$ across the $T$ rounds. The time needed to inflate and
deflate are dominated by $O(Tkn)$, establishing Lemma 7.

##  Appendix H Locating $\lambda$ and Describing the Optimal Adversarial
Strategy

###  H.1. Proof of Theorem 8

We first begin with a toy version that says if the rewards from
$\operatorname{LinearCSSPA}(\alpha,\beta,\lambda_{1},T,k)$ and
$\operatorname{LinearCSSPA}(\alpha,\beta,\lambda_{2},T,k)$ are close, then
$\lambda_{1}-\lambda_{2}$ is small.

###### Lemma 0.

Suppose that $\lambda_{1}$ and $\lambda_{2}$ are such that

$|\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{1}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{1}))-\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{2}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{2}))|\leq\zeta\text{.}$

Then, $|\lambda_{1}-\lambda_{2}|\leq\zeta$.

###### Proof.

Without loss of generality, assume that $\lambda_{1}\leq\lambda_{2}$. Assume
for contradiction that $\lambda_{2}-\lambda_{1}>\zeta$.

As a thought experiment, let the adversary ape
$\pi^{\operatorname{OPT}}_{T,k}(\lambda_{2})$ in
$\operatorname{LinearCSSPA}(\alpha,\beta,\lambda_{1},T,k)$. Since the rewards
are linear, the adversary earns
$\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{2}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{2}))$
plus the savings in entry fee from paying just $\lambda_{1}<\lambda_{2}$.
Since the adversary participates in at least one round, the adversary saves
more than $(\lambda_{2}-\lambda_{1})$ in entry fee, and thus, the total reward
is at least
$\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{2}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{2}))+(\lambda_{2}-\lambda_{1})$.
However,
$\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{1}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{1}))-\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{2}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{2}))\leq\zeta$
and $(\lambda_{2}-\lambda_{1})\geq\zeta$. The adversary makes strictly more
than
$\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{1}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{1}))$
in $\operatorname{LinearCSSPA}(\alpha,\beta,\lambda_{1},T,k)$. This is a
contradiction, since $\pi^{\operatorname{OPT}}_{T,k}(\lambda_{1})$ is the
optimal strategy. ∎

In Theorem 8, we show that even if we estimate
$\operatorname{LinearCSSPA}(\alpha,\beta,\lambda_{1},T,k)$ up to an error of
$\delta$ then, the true reward in
$\operatorname{LinearCSSPA}(\alpha,\beta,\lambda_{2},T,k)$ and our estimated
rewards can still not be close unless $\lambda_{1}-\lambda_{2}$ is small.

###### Proof of Theorem 8.

From Theorem 5, with probability at least
$1-\Big{(}2\,T\,\gamma+T\,\tfrac{e^{-\omega
n}}{\omega}\,\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}\Big{)}$,

$\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{1}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{1}))\in\Big{[}\mathbb{E}_{s\sim\overline{\operatorname{D}}_{T,k}(\lambda_{1})}[s],\mathbb{E}_{s\sim\underline{\operatorname{D}}_{T,k}(\lambda_{1})}[s]\Big{]}$

Thus,
$|\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{1}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{1}))-r|\leq\delta$
and
$|r-\operatorname{Rew}^{\operatorname{Lin}}_{\lambda_{2}}(\pi^{\operatorname{OPT}}_{T,k}(\lambda_{2}))|\leq\zeta$.
The claim follows from the triangle inequality and Lemma 1. ∎

###  H.2. The Optimal Adversarial Strategy

Once we (approximately) locate the optimal reward $\lambda^{*}$ of the
$k$-scored adversary in $\operatorname{CSSPA}(\alpha,\beta)$ that terminates
in $T$ rounds, the description of a near optimal adversarial strategy becomes
succinct. The adversary pretends that it is $k$-scored and is participating in
$\operatorname{LinearCSSPA}(\alpha,\beta,\lambda,T)$ and plays the actions
recommended by $\pi^{\operatorname{OPT}}_{T,k}(\lambda^{*})$. In particular,
the adversary resets $\operatorname{CSSPA}(\alpha,\beta)$ with frequency at
least $T$ and behaves like a $k$-scored adversary. The adversary gets a zero
reward in $\operatorname{LinearCSSPA}(\alpha,\beta,\lambda,T)$, which
corresponds to a reward $\lambda^{*}$ in $\operatorname{CSSPA}(\alpha,\beta)$.

The reward from the above strategy approaches the optimal reward as the error
in locating $\lambda^{*}$ reduces. The loss in the reward from being
$k$-scored and from resetting $\operatorname{CSSPA}(\alpha,\beta)$ every $T$
rounds would approach zero as $T\xrightarrow{}\infty$ and
$k\xrightarrow{}\infty$.

##  Appendix I Summary

###  I.1. A Summary of the Simulation

$\text{FiniteSampleAddLayer}(\alpha,\beta,\lambda,k,t,n,\pi^{\operatorname{OPT}}_{T,k},\gamma,\omega,\epsilon,\eta,\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k})$:
(1)
$\text{Precompute}(\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k},\epsilon,\eta,t)$:
(a) Construct the pdf $d$ of
$\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}$ up to a discretization
error $\epsilon$. (b) Construct the cdf
$\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}$ recursively using the
following recursion up to a discretization error $\epsilon$.
$\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}(\theta-\epsilon)=Pr(r_{0}\leq\theta-\epsilon)=\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}(\theta)-d(\theta)$
(c) Construct the right tail of the expectation
$E(\theta)=\mathbb{E}_{r_{0}\sim D}[r_{0}\times\mathbbm{1}(r_{0}>\theta)]$
recursively by $E(\theta-\epsilon)=E(\theta)+(\theta-\epsilon)\times
d(\theta-\epsilon)$ (d) Compute
$E_{\max}(\theta)=\theta\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k}+E(\theta)$
for all $-t\lambda\leq\theta\leq t(1-\lambda)$ up to a discretization error
$\epsilon$ (e) If $\beta\neq 1,0$: For $-t\lambda\leq
g(i^{*},\vec{c}_{-0},\vec{r}_{-0})\leq t(1-\lambda)$ in steps of $\epsilon$:
(i) For $0\leq c\leq 1$ in steps of $\eta$ •
$G(c,g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))=\eta\times\sum_{\zeta=0,\text{ in
steps of
}\eta}^{c}\zeta^{\frac{1-\beta}{\beta}}E_{\max}(\frac{g(i^{*},\vec{c}_{-0},\vec{r}_{-0})}{\zeta^{\frac{1-\beta}{\beta}}})$
(2) For $1\leq\ell\leq n$. (a)
$\text{DrawAdv}(\alpha,\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k})$:
• Sample $\vec{r}_{-0}$: Draw $k$ rewards $r_{1},r_{2},\dots,r_{k}$ i.i.d from
$\operatorname{D}^{\operatorname{OPT}}_{t-1}$. • Sample $\vec{c}_{-0}$: Draw
$k$ scores $c_{1},c_{2},\dots,c_{k}$ of adversarial wallet as follows. Draw
$c_{1}\xleftarrow{}\operatorname{exp}(\alpha)$ and
$c_{i+1}\xleftarrow{}c_{i}+\operatorname{exp}(\alpha)$ (a fresh sample for
each $i$) for $1\leq i\leq k-1$. For convenience, set $c_{k+1}=\infty$. •
Return $(\vec{r}_{-0},\vec{c}_{-0})$. (b)
$\text{SampleFromPrecompute}(\alpha,\beta,\lambda,\vec{c}_{-0},\vec{r}_{-0},\pi^{\operatorname{OPT}}_{T,k},G,\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k})$:
(i) For $1\leq i^{*}\leq k$: • Compute
$g(i^{*},\vec{c}_{-0},\vec{r}_{-0})=\max_{1\leq i\leq
i^{*}}{e^{-c_{i}\cdot(1-\beta)(1-\alpha)}}(1+r_{i})$ (ii) If $\beta=0$: •
return $s_{\ell}=g(k,\vec{c}_{-0},\vec{r}_{-0})-\lambda$ (iii) Else if
$\beta=1$: • For $1\leq i^{*}\leq k$:
$f(i^{*},\vec{c}_{-0},\vec{r}_{-0})=E_{\max}(g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))\big{[}e^{-c_{i^{*}}\,(1-\alpha)}-e^{-c_{i^{*}+1}\,(1-\alpha)}\big{]}$
• return
$s_{\ell}=\sum_{i^{*}=1}^{k}f(i^{*},\vec{c}_{-0},\vec{r}_{-0})-\lambda$ (iv)
$\beta\neq 0,1$: • For $1\leq i^{*}\leq k$:
$f(i^{*},\vec{c}_{-0},\vec{r}_{-0})=G(c_{i^{*}+1},g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))-G(c_{i^{*}},g(i^{*},\vec{c}_{-0},\vec{r}_{-0}))$
• return
$s_{\ell}=\sum_{i^{*}=1}^{k}f(i^{*},\vec{c}_{-0},\vec{r}_{-0})-\lambda$ (3)
$\tilde{\operatorname{D}}^{\operatorname{OPT}}_{t,k}$ be the uniform
distribution over $\big{(}s_{\ell}\big{)}_{1\leq\ell\leq n}$ (in descending
order) (4) Inflate while computing the upper bound and deflate while computing
the lower bound. •
$\text{Deflate}(n,\gamma,\tilde{\operatorname{D}}^{\operatorname{OPT}}_{t,k},t):$
(a) Delete the largest $n\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$ samples from
$\tilde{\operatorname{D}}^{\operatorname{OPT}}_{t,k}$ (b) Append
$n\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$ copies of $-\lambda$ to
$\tilde{\operatorname{D}}^{\operatorname{OPT}}_{t,k}$ •
$\text{Inflate}(n,\gamma,\omega,\tilde{\operatorname{D}}^{\operatorname{OPT}}_{t,k},t)$:
(a) Delete the smallest $n\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$ samples from
$\tilde{\operatorname{D}}^{\operatorname{OPT}}_{t,k}$ (b) Append $\omega n$
copies of $t\,(1-\lambda)$ to
$\tilde{\operatorname{D}}^{\operatorname{OPT}}_{t,k}$ (c) For
$1\leq\ell<\frac{n}{\omega\,n}\cdot\sqrt{\tfrac{\ln\gamma^{-1}}{2n}}$: –
Append $\omega n$ copies of $s_{\ell}$ (5) Return
$\hat{\operatorname{D}}^{\operatorname{OPT}}_{t,k}$ to be the uniform
distribution over $\\{s_{1},s_{2},\dots,s_{n_{t}}\\}$.
$\text{TruncatedSimulate}(\alpha,\beta,\lambda,T,k,\pi^{\operatorname{OPT}}_{T,k},\gamma,\omega,\epsilon,\eta)$:
(1) Initialize $\hat{\operatorname{D}}^{\operatorname{OPT}}_{0,k}$ to be the
point-mass distribution at $0$. (2) For $1\leq t\leq T$: (a)
$\hat{\operatorname{D}}^{\operatorname{OPT}}_{t,k}=\text{FiniteSampleAddLayer}(\alpha,\beta,\lambda,k,t,n,\pi^{\operatorname{OPT}}_{T,k},\gamma,\omega,,\epsilon,\eta,\hat{\operatorname{D}}^{\operatorname{OPT}}_{t-1,k})$.
(3) Return
$\mathbb{E}_{s\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{T,k}}[s]$.

###  I.2. A Summary of Notations and Functions in the Simulation

Notation | Description
---|---
$\alpha$ | The fraction of stake held by the adversary
$\beta$ |  The fraction of honest stake held by $B$
$\lambda$ |  The per-round entry fee the adversary has to pay to participate in $\operatorname{LinearCSSPA}(\alpha,\beta,\lambda)$
$T$ |  We simulate $\operatorname{LinearCSSPA}(\alpha,\beta,\lambda,T)$
$k$ |  We simulate a $k$-scored adversary
$\pi^{\operatorname{OPT}}_{T,k}$ |  The optimal adversarial strategy in $\operatorname{LinearCSSPA}(\alpha,\beta,\lambda,T)$ for a $k$-scored adversary
$\hat{\operatorname{D}}^{\operatorname{OPT}}_{t,k}$ |  The estimated distribution of optimal rewards in $\operatorname{LinearCSSPA}(\alpha,\beta,\lambda,t)$ for a $k$-scored adversary
$\gamma$ |  Probability of estimation error from the DKW inequality in inflate/deflate
$\omega$ |  A small quantile gets duplicated $\omega n$ times while inflating
$\epsilon$ |  Discretization parameter of the reward distributions
$\eta$ |  Discretization parameter to compute Riemann sums

Table 1. A summary of notations

Function | Description
---|---
$\text{AddLayer}(\alpha,\beta,\lambda,\pi,D)$ |  Given a distribution $D$ of rewards achieved by playing $\pi$ in the last $t-1$ rounds, $\text{AddLayer}(\alpha,\beta,\lambda,\pi,D)$ computes the distribution of rewards won by playing $\pi$ in the last $t$ rounds of $\operatorname{LinearCSSPA}(\alpha,\beta,\lambda,t)$
$\text{DrawAdv}(\alpha,D)$ |  Given the reward distribution $D$, $\text{DrawAdv}(\alpha,D)$ samples rewards and scores for adversarial wallets
$\text{sample}(\alpha,\beta,\lambda,\vec{c}_{-0},\vec{r}_{-0},\pi,D)$ |  Conditional on the rewards and scores of the adversary’s wallets, $\text{sample}(\alpha,\beta,\lambda,\vec{c}_{-0},\vec{r}_{-0},\pi,D)$ computes the reward of the adversary, in expectation over $B$ and $C$’s scores and the reward from letting $B$ win.
$\text{Inflate}(n,\gamma,\omega,\operatorname{D},t)$ |  For an input distribution $\operatorname{D}$, $\text{Inflate}(n,\gamma,\operatorname{D},t)$ deletes the smallest samples and replaces them with samples corresponding to a small quantile
$\text{Deflate}(n,\gamma,\operatorname{D},t)$ |  For an input distribution $\operatorname{D}$, $\text{Deflate}(n,\gamma,\operatorname{D},t)$ deletes the largest samples and replaces them with the infimum of the distribution
$\text{Precompute}(\operatorname{D},\epsilon,\eta,t)$ |  For an input distribution $\operatorname{D}$, $\text{Precompute}(D,\epsilon,\eta,t)$ sets up the pre-computations required to speed up the sampling procedure
$\text{Simulate}(\alpha,\beta,\lambda,\pi^{\operatorname{OPT}}_{T,k})$ |  Executes $T$ iterations of $\text{AddLayer}(\alpha,\beta,\lambda,\pi^{\operatorname{OPT}}_{T,k})$ and returns $\mathbb{E}_{s\sim\hat{\operatorname{D}}^{\operatorname{OPT}}_{T,k}}[s]$

Table 2. A summary of functions used in the simulation. Functions with similar
names across different variants of the simulation have similar functions. For
example,
$\text{AddLayer}(\alpha,\beta,\lambda,\pi^{\operatorname{OPT}}_{T,k},D)$ and
$\text{FiniteSampleAddLayer}(\alpha,\beta,\lambda,k,t,n,\pi^{\operatorname{OPT}}_{T,k},\gamma,\omega,D)$
have the sample functionality (adding an extra layer in the simulation), but
in different variants.


# A reduced variational approach for searching cycles in high-dimensional
systems

Ding Wang1 Yueheng Lan1,2<EMAIL_ADDRESS>1School of Science, Beijing
University of Posts and Telecommunications, Beijing 100876, China.
2State Key Lab of Information Photonics and Optical Communications,Beijing
University of Posts and Telecommunications, Beijing 100876, China

###### Abstract

Searching recurrent patterns in complex systems with high-dimensional phase
spaces is an important task in diverse fields. In the current work, an
improved scheme is proposed to accelerate the recently designed variational
approach for finding periodic orbits in systems with chaotic dynamics based on
the existence of inertial manifold widely observed in various spatially
extended systems, especially those with high dimensions. On the premise of
keeping exponential convergence of the variational method, an effective loop
evolution equation is derived to greatly reduce the storage and computing
time. With repeated modification of local coordinates and evolution of the
guess loop being carried out alternately, the rapid convergence and the
stability of the reduction scheme are effectively achieved. The dimension of
local coordiante subspaces is generally larger than the number of nonnegative
Lyapunov exponents to ensure the exponential convergence. The proposed scheme
is successfully demonstrated on several well-known examples and expected to
supply a powerful tool in the exploration of high-dimensional nonlinear
systems.

reduced variational method, the complex system

##  I Introduction

Complex phenomena in the fields of geophysics, space physics and fluid
mechanics are usually described by continuous-time dynamical systems, in which
infinitely many unstable periodic orbits (cycles) play important roles in
their characterization and analysis. In the study of turbulence[1, 2],
unstable periodic orbits (UPOs) extracted from fluid flows can well describe
the characteristics of turbulence, such as mean velocity profile, the time
development of coherent structures, intermittency or other statistical
properties. In a plethora of contexts, periodic orbits could be conveniently
used to study the long-term behavior of a nonlinear system [3], for which
cycle expansions are very efficient in computing average values of physical
observables [4, 5, 6] on the strange attractors, where cycles are ordered
hierarchically in terms of their topological length or stability. Therefore,
it is of great significance to develop efficient numerical methods and
techniques for locating periodic orbits.

In literature, various numerical algorithms of identifying periodic orbits
have been proposed so far[7, 8, 9, 10, 11, 12, 13, 14], among which the best
known is probably the Newton-Raphson method. The method and its variants[8, 9]
show great convenience and efficiency in finding short cycles in a large
variety of dynamical systems, which however for UPOs with longer periods
requires good initial guesses. Even still, it is likely to fail in finding
unstable periodic orbits in high-dimensional systems[15], because the
cumulative error in the long-time evolution increases exponentially with the
evolution time. The multiple shooting algorithm is able to overcome these
difficulties but a set of Poincar$\acute{\text{e}}$ sections need to be
introduced to divide the whole trajectory into short segments, so as to
contain the error of the long-time integration. Nevertheless, it seems hard to
select suitable Poincar$\acute{\text{e}}$ sections in high dimensions, which
should intersect all the neighboring orbits transversely. To remove this
trouble, the damped Newton-Raphson-Mees algorithm is proposed in Ref.[14], but
the types of UPOs that can be detected by the method are limited. An
alternative scheme is proposed in Ref.[13], which describes a periodic orbit
with a loop of discrete points. A variational equation derived from the
original dynamics will drive the guess loop to a genuine periodic orbit
exponentially. Because of the topological constraint embodied in the loop
representation, irrespective of the original dynamics, the method shows
excellent computational stability.

However, the variational method requires the storage and inversion of an
$(Nd+1)\times(Nd+1)$ matrix ($d$ is the dimension of system and $N$ is the
number of lattice points on the guess loop), which is a huge burden when
exploring complex systems. An automatic mesh allocation scheme[16] and
additional techniques on storage optimization are designed to alleviate these
troubles to some extent in search of connecting orbits[17] or periodic orbits
near singularity[18]. However, these strategies can not fundamentally break
the bottleneck of the computational complexity in high dimensions.

In the current work, to extend application of the variational method to more
stringent situation, we propose a novel scheme to greatly reduce the
computation load based on the consideration below. If a cycle is stable, any
point in its neighborhood will evolve to it following the equations of motion.
Therefore, for unstable cycles, only the unstable or neutral directions need
to be controlled during evolution. The good news is that in most high-
dimensional systems, the number of expanding directions is much smaller than
the phase space dimension. Hence, it is possible to build a numerical scheme
which cares only about these unstable directions so as to accelerate the
computation. The framework of the variational algorithm is especially good for
this purpose. Along the guess loop, a fictitious evolution could be designed
(see below) to dig out all important directions, the number of which is, say,
$n$. With $n\ll d$, the complexity of the matrix computation and hence of the
whole scheme is much reduced!

The paper is organized as follows. In Sec. II, The variational principle is
briefly introduced. In Sec. III.2, The detailed derivation of the reduced
variational equation and the construction and evolution of local coordinates
are given. In Sec. IV, three examples are used to demonstrate the validity and
application details of the reduced variational approach. Our computations are
summarized and discussions are made in the final Section V.

##  II the variational principle of cycle searching

For a general $d$-dimensional dynamical system, its time evolution could be
depicted by a set of ordinary differential equations (ODEs):

$\dot{\bm{x}}(t)=\bm{v}(\bm{x}(t)),$ (1)

where $\bm{x}\in\mathbb{R}^{d}$, $t\in\mathbb{R}$ and
$\bm{v}(\bm{x}):\mathbb{R}^{d}\to\mathbb{R}^{d}$ is a smooth function defined
in the phase space. A periodic orbit ( or interchangeably cycle) with period
$T$ is a trajectory $\bm{x}(t)$ of Eq. (1) that satisfies the condition

$f^{T}(\bm{x})=\bm{x},$ (2)

where $f^{T}(\bm{x})$ gives the new position of $\bm{x}$ under the system
dynamics Eq. (1) after a time interval $T>0$. The variational method and its
variants have been widely used to locate cycles in many systems[13, 18, 17]
and show impressive robustness and convergence, which deploy a loop in the
phase space being driven to a true cycle by the fictitious evolution equation

$\frac{\partial^{2}\tilde{\bm{x}}}{\partial
s\partial\tau}-\lambda\frac{\partial\bm{v}}{\partial\bm{x}}\frac{\partial\tilde{\bm{x}}}{\partial\tau}-\bm{v}\frac{\partial\lambda}{\partial\tau}=\lambda\bm{v}-\tilde{\bm{v}},$
(3)

where $\tilde{\bm{x}}$ marks the representative points of the loop and is
parameterized by $s~{}(s\in[0,2\pi])$,
$\tilde{\bm{v}}=\frac{\partial\tilde{\bm{x}}}{\partial
s}~{}(\tilde{\bm{x}}\in\mathbb{R}^{d})$ is the loop velocity, $\bm{v}$ is the
velocity field of the dynamical system given by Eq. (1) and $\tau$ is the
fictitious time which records the evolution from an initial guess loop to the
desired periodic orbit. The parameter $\lambda$ is independent of $s$ but is a
function of $\tau$ which is used to match the magnitude of $\bm{v}$ and
$\tilde{\bm{v}}$, or in other words, to adjust the period. We may rewrite Eq.
(3) as

$\frac{\partial}{\partial\tau}(\tilde{\bm{v}}-\lambda\bm{v})=\lambda\bm{v}-\tilde{\bm{v}},$
(4)

the solution of which could be formally written as

$\tilde{\bm{v}}-\lambda\bm{v}=e^{-\tau}(\tilde{\bm{v}}-\lambda\bm{v})|_{\tau=0},$
(5)

showing that the error $\tilde{\bm{v}}-\lambda\bm{v}$ exponentially decreases
with $\tau$ and the loop exponentially approaches the periodic orbit.

In a discretization of the loop, the vector $\tilde{\bm{v}}$ is calculated by
a five-point approximation scheme for numerical stability:

$\displaystyle\tilde{\bm{v}}(s_{i})$
$\displaystyle\equiv\frac{1}{12h}[8\tilde{\bm{x}}(s_{i+1})-8\tilde{\bm{x}}(s_{i-1})+\tilde{\bm{x}}(s_{i-2})-\tilde{\bm{x}}(s_{i+2})]$
$\displaystyle=\hat{\bm{D}}\tilde{\bm{x}}(s_{i})$ (6)

in which $\hat{\bm{D}}$ is a matrix operator with the periodic condition

$\hat{\bm{D}}=\frac{1}{12h}\begin{pmatrix}~{}~{}0&~{}~{}8&-1&~{}&~{}&~{}&~{}~{}1&-8\\\
-8&~{}~{}0&~{}~{}8&-1&~{}&~{}&~{}&~{}~{}1\\\
~{}~{}1&-8&~{}~{}0&~{}~{}8&-1&~{}&~{}&~{}\\\
~{}&~{}~{}1&-8&~{}~{}0&~{}~{}8&-1&~{}&~{}\\\
~{}&~{}&\ddots&\ddots&\ddots&\ddots&\ddots&~{}\\\
~{}&~{}&~{}&~{}~{}1&-8&~{}~{}0&~{}~{}8&-1\\\
-1&~{}&~{}&~{}&~{}~{}1&-8&~{}~{}0&~{}~{}8\\\
~{}~{}8&-1&~{}&~{}&~{}&~{}~{}1&-8&~{}~{}0\end{pmatrix}$ (7)

where every element is a $d\times d$ diagonal matrix and $h=2\pi/N$, $N$ being
the number of lattice points. The discretized version of Eq. (3) with a
fictitious time interval $\delta\tau$ is

$\begin{pmatrix}\hat{\bm{A}}&~{}-\hat{\bm{v}}\\\
\hat{\bm{a}}&~{}~{}~{}0\end{pmatrix}\begin{pmatrix}\delta\hat{\bm{x}}\\\
\delta\lambda\end{pmatrix}=\delta\tau\begin{pmatrix}\lambda\hat{\bm{v}}-\hat{\tilde{\bm{v}}}\\\
0\end{pmatrix}$ (8)

with

$\hat{\bm{A}}=\hat{\bm{D}}-\lambda\text{diag}[A_{1},~{}A_{2},...,A_{N}]$

where $A_{n}$ is the velocity gradient matrix $A_{n}=\frac{\partial
v_{i}}{\partial v_{j}}|_{s=s_{n}}$,
$\hat{\bm{v}}=(\bm{v}_{1},\bm{v}_{2},...,\bm{v}_{n})^{t}$ and
$\hat{\tilde{\bm{v}}}=(\tilde{\bm{v}}_{1},\tilde{\bm{v}}_{2},...,\tilde{\bm{v}}_{n})^{t}$
are the velocity vector of the original flow and the loop velocity,
respectively. $\hat{\bm{a}}$ is an $Nd$-dimensional row vector which imposes a
gauge-fixing condition on the coordinate corrections $\delta\hat{\bm{x}}$ in
order to remove the neutral direction along which the cycle is invariant upon
shifting all the points. In each iteration, we set up and numerically solve
Eq. (8) and update the coordinates of the guess loop which is supposed to
approach exponentially a periodic orbit if it exists.

##  III A reduced numerical Scheme

Most computation load originates from the matrix
$\tiny{\begin{pmatrix}\hat{\bm{A}}&~{}-\hat{\bm{v}}\\\
\hat{\bm{a}}&~{}~{}~{}0\end{pmatrix}}$ in Eq. (8) which is
$(Nd+1)\times(Nd+1)$, involving the dimension $d$ of the phase space and the
number of lattice points $N$. With the increase of $d$ and $N$, the required
storage and computing time quickly increases, which is a major bottleneck
preventing a wide application of the variational method in high-dimensional
systems. In Ref. [18], an automatic allocation scheme of lattice points is
adopted to minimize $N$. However, in a complex system, the dimension $d$ of
the system is much greater than $N$. Only reducing the number of lattice
points is far from enough, and thus we try to accelerate the variational
method by reducing the effective dimension of the local coordinate system.

###  III.1 Reduction by projection

In practice, not all directions are equally important for orbit adjustment
when approaching a periodic orbit from a guess loop. On a hyperplane
perpendicular to the periodic orbit, the vicinity of an orbit point is
stretched in a few directions and compressed in others when moving along the
orbit. The correction in the compressed directions is automatically made
during the orbit evolution and the deviation in the stretching or neutral
directions has to be corrected by the variational scheme. Therefore, it is
essential to find these important directions and rewrite the original
variational equation in reduced coordinate frames.

Figure 1:  The local coordinate frames $\\{\bm{e}_{k}(s_{i})\\}$
$(k=1,2,...,n;~{}i=1,2,...,N$ and $n\leq d)$. along the orbit parametrized by
$s_{i}$, where $S_{i}$ represents the hyperplane spanned by the vectors
$\bm{e_{k}}(s_{i})~{}(k\leq d)$.

Assuming that along the important directions, we already built a family of
local coordinate frames $\\{\bm{e}_{k}(s_{i})\\}$
$(k=1,2,...,n;~{}i=1,2,...,N$ and $n\leq d)$ for the projection ( Fig.1), the
reduced variational equation may be derived as follows.

Firstly, we write an approximation of Eq.(II) in the subspace
$\\{\bm{e}_{k}\\}$,

$\displaystyle\tilde{v}_{j}(s_{i})$
$\displaystyle=\frac{1}{12h}[8\tilde{x}_{j}(s_{i+1})\bm{e}_{j}(s_{i+1})-8\tilde{x}_{j}(s_{i-1})\bm{e}_{j}(s_{i-1})+\tilde{x}_{j}(s_{i-2})\bm{e}_{j}(s_{i-2})-\tilde{x}_{j}(s_{i+2})\bm{e}_{j}(s_{i+2})],$
(9)

and a dot-product of both sides of Eq.(3) with $\bm{e}_{k}(s_{i})$ results in

$\bm{e}_{k}(s_{i})\cdot\bigg{[}\frac{\partial^{2}\tilde{\bm{x}}}{\partial
s\partial\tau}-\lambda\frac{\partial\bm{v}}{\partial\bm{x}}\frac{\partial\tilde{\bm{x}}}{\partial\tau}-\bm{v}\frac{\partial\lambda}{\partial\tau}\bigg{]}=\bm{e}_{k}(s_{i})\cdot\bigg{[}\lambda\bm{v}-\tilde{\bm{v}}\bigg{]}.$
(10)

By feeding Eq. (9) into Eq. (10), we obtain the projection of the first term
in Eq.(10),

$\bm{e}_{k}(s_{i})\cdot\Delta\tilde{\bm{v}}(s_{i})=\frac{1}{12h}\bigg{[}\Delta\tilde{x}_{i-2}^{j}\cdot(\bm{e}^{k}_{i}\cdot\bm{e}^{j}_{i-2})-8\Delta\tilde{x}_{i-1}^{j}(\bm{e}^{k}_{i}\cdot\bm{e}^{j}_{i-1})+8\Delta\tilde{x}_{i+1}^{j}\cdot(\bm{e}^{k}_{i}\cdot\bm{e}^{j}_{i+1})-\Delta\tilde{x}_{i+2}^{j}\cdot(\bm{e}^{k}_{i}\cdot\bm{e}^{j}_{i+2})\bigg{]},$
(11)

where $\Delta\tilde{x}^{j}_{i}$ and $\bm{e}^{j}_{i}$ labels
$\Delta\tilde{x}_{j}(s_{i})$ and $\bm{e}_{j}(s_{i})$ respectively. According
to Eq. (11), every constant block in the differential operator $\hat{\bm{D}}$
in Eq. (7) is replaced by an $n\times n$ matrix which is defined by
$\bm{e}^{k}_{i}\cdot\bm{e}^{j}_{i^{\prime}}$, where $(k,j)$ marks the block
position and $(i,i^{\prime})$ the position within the $(k,j)-$block. The
second term
$\frac{\partial\bm{v}}{\partial\bm{x}}\frac{\partial\tilde{\bm{x}}}{\partial\tau}$
in Eq.(10) can be written as

$\displaystyle\bm{e}_{k}(s_{i})\cdot\Delta\bm{v}(s_{i})$
$\displaystyle=\bm{e}_{k}(s_{i})\cdot\bigg{[}\bm{v}\bigg{(}\tilde{\bm{x}}_{i}+\sum\limits_{j}\Delta
x_{ij}\bm{e}_{j}\small(s_{i}\small)\bigg{)}-\bm{v}(\tilde{\bm{x}}_{i})\bigg{]}$
$\displaystyle=\bm{e}_{k}(s_{i})\cdot\sum\limits_{p}\bigg{[}\hat{\bm{e}}_{p}v_{p}\big{(}\tilde{\bm{x}}_{i}+\sum\limits_{j}\Delta
x_{ij}\bm{e}_{j}\small(s_{i}\small)\big{)}-\hat{\bm{e}}_{p}v_{p}(\tilde{\bm{x}}_{i})\bigg{]}$
$\displaystyle=\bm{e}_{k}(s_{i})\cdot\sum\limits_{p}\bigg{[}\hat{\bm{e}}_{p}v_{p}(\tilde{\bm{x}}_{i})+\hat{\bm{e}}_{p}\sum\limits_{j}\big{(}\frac{\partial
v_{p}}{\partial x_{j}}\Delta
x_{ij}\big{)}-\hat{\bm{e}}_{p}v_{p}(\tilde{\bm{x}}_{i})\bigg{]}$
$\displaystyle=\sum\limits_{p}\bigg{[}{e}_{kp}(s_{i})\sum\limits_{j}\big{(}\frac{\partial
v_{p}}{\partial x_{j}}\Delta x_{ij}\big{)}\bigg{]}$
$\displaystyle=\sum\limits_{j}\frac{\partial v_{k}}{\partial x_{j}}\Delta
x_{ij}$ (12)

where $k\,,j$ marks the new coordinates and $\\{\hat{\bm{e}}_{p}\\}$ make up
the original Cartesian coordinate system which are constant at every point in
the phase space. $v_{k}=\bm{\bm{v}}\cdot\bm{e}_{k}(s_{i})$ ($\\{\bm{e}_{k}\\}$
is an orthogonal basis),
$e_{kp}(s_{i})=\bm{e}_{k}(s_{i})\cdot\hat{\bm{e}}_{p}$, and
$\Delta\tilde{\bm{x}}(s_{i})=\sum\limits_{j}\Delta
x_{ij}\bm{e}_{j}(s_{i})=\Delta\tilde{\bm{x}}_{i}$. Therefore, the size of the
velocity gradient matrix $\frac{\partial\bm{v}_{k}}{\partial\bm{x}_{j}}$,
which is obtained by numerical differentiation, is much smaller than that of
the original velocity gradient $A=\partial v/\partial x$ of Eq.(3) if $n\ll
d$. About other terms in Eq.(10), the vectors can be just projected into the
local subspace $\\{\bm{e}_{k}\\}$ and the matrix
$\tiny{\begin{pmatrix}\vec{\bm{A}}&~{}-\vec{\bm{v}}\\\
\vec{\bm{a}}&~{}~{}~{}0\end{pmatrix}}$ in Eq. (8) is then reduced to
$(Nn+1)\times(Nn+1)$. Because $n\ll d$ actually holds in many nonlinear
systems, the reduced matrix order $Nn+1$ will be much smaller than the
original one $Nd+1$. Therefore the storage and computing time required to deal
with the reduced form of Eq. (8) are greatly cut down. The larger the
dimension $d$, the more prominent the benefits of this reduction.

###  III.2 Construction and evolution of the local coordinates

In this section, we will discuss how to build the local coordinate system
$\\{\bm{e}_{k}\\}$ mentioned above at each lattice point $s_{i}$. However,
before that, Let’s first explain the concept of pseudo-evolution. Two adjacent
points on the loop are approximately related by the evolution
$\bm{x}(s_{i})\longrightarrow\bm{x}(s_{i+1})$ with Eq. (1) if the guess is
reasonably good. Similarly, the corresponding nearby points are assumed
evolving as
$\bm{x}(s_{i})+\bm{h}(s_{i})\longrightarrow\bm{x}(s_{i+1})+\bm{h}(s_{i+1})$
governed by Eq.(1) either, where $\bm{h}(s_{i})$ and $\bm{h}(s_{i+1})$
represent respectively tiny offset vectors at $\bm{x}(s_{i})$ and
$\bm{x}(s_{i+1})$. Then we have approximately

$\bm{h}(s_{i+1})\longrightarrow\bm{h}(s_{i})+[\bm{v}(\bm{x}_{i}+\bm{h}(s_{i}))-\bm{v}(\bm{x}_{i})]\Delta
t,$ (13)

where $\Delta t$ is the time interval between points $s_{i}$ and $s_{i+1}$.
Eq. (13) is the pseudo-evolution equation of points near the guess loop, which
obviously is rough when the guess loop is far away from the periodic orbit we
are looking for. As the loop gradually approaches the periodic orbit, it
becomes more and more accurate.

In order to decently depict the evolution along the orbit, we employ a second-
order numerical scheme based on Eq. (13),

$\displaystyle\bm{v}_{1}=[\bm{v}(\bm{x}_{i})+\bm{v}(\bm{x}_{i+1})]/2,$
$\displaystyle\bm{v}_{2}=[\bm{v}(\bm{x}_{i}+\bm{h}_{j}(s_{i}))+\bm{v}(\bm{x}_{i+1}+\bm{h}_{j}(s_{i+1}))]/2,$
(14)
$\displaystyle\bm{h}^{\prime}_{j}(s_{i+1})\longrightarrow\bm{h}_{j}(s_{i})+(\bm{v}_{2}-\bm{v}_{1})\Delta
t,$

where $\bm{h}_{j}(s_{i})$ is a tiny vector in the $j$-$th$ direction
$\bm{e}_{j}(s_{i})$, $\bm{h}^{\prime}_{j}(s_{i+1})$ is the tiny vector at
$s_{i+1}$ after one step of evolution, while $\bm{h}_{j}(s_{i+1})$ is the
vector at point $s_{i+1}$ before the evolution. In order to ensure that the
new coordinate system is orthogonal, after each evolution it is necessary to
carry out the Schmidt orthogonalization and normalization of
$\\{\bm{h}_{j}^{\prime}(s_{i})\\}$ to get the new set
$\\{\bm{e}_{j}(s_{i})\\}$ at each $s_{i}$,

$\displaystyle\bm{h}_{j}^{\prime}(s_{i})\longrightarrow\bm{h}_{j}^{\prime}(s_{i})-\sum_{k=1}^{j-1}c_{k}\bm{e}_{k}(s_{i}),$
$\displaystyle\bm{e}_{j}(s_{i})\longrightarrow\bm{h}_{j}^{\prime}(s_{i})/|\bm{h}_{j}^{\prime}(s_{i})|,$

where $c_{k}=\bm{e}_{k}(s_{i})\cdot\bm{h}^{\prime}_{j}(s_{i})$. As shown in
Fig. 2, the evolution of the local coordinate axis in the $j$-$th$ direction
is plotted. For each evolution, a very small vector $\bm{h}_{j}$ (black arrow)
in some $j$-th direction at each lattice point $s_{i}$ evolves to
$\bm{h}_{j}^{\prime}(s_{i+1})$ (red arrow) at $s_{i+1}$. After orthomalizing
the vector $\bm{h}_{j}^{\prime}(s_{i+1})$, we get a new direction vector
$\bm{e}_{j}(s_{i+1})$ (not displayed) at $s_{i+1}$. Repeating the procedure in
all involved directions we will finally build up a new coordinate frame at
$s_{i+1}$.

Figure 2:  The evolution of local coordinate frame. The thick black solid line
represents the guess periodic orbit $\bm{X}(s_{i})$, and the dotted line is an
adjacent curve $\bm{X}_{near}(s_{i})=\bm{X}(s_{i})+\bm{h}_{j}$ near the
periodic orbit, and $|\bm{h}_{j}|=10^{-3}$. $s_{i}$ is the curvilinear
coordinate of the point on the guess orbit. $\bm{h}^{\prime}_{j}(s_{i+1})$
(red arrow) at $s_{i+1}$ evolves from $\bm{h}_{j}$ (black arrow) at $s_{i}$.

What is a proper dimension of the local coordinate system? Obviously,
according to previous arguments it should not be less than the number of
positive Lyapunov exponents of the system. The Lyapunov exponent equals $0$ in
the direction of the velocity field, but the sampling points $x(s_{i})$ will
approach the desired periodic orbit only along stable directions. Therefore,
we include the velocity vector in the reduced coordinate subspace, as well as
$2\sim 3$ other not-so-unstable or neutral directions to accelerate the
convergence. To start with, at each lattice point we may select the same first
few coordinate directions in the full phase space to build the initial local
frame. Obviously, such a coordinate system is not in the stretching direction
of the orbit, and we should have evolved it for some time to achieve that.
Nevertheless, it does not hinder us from computing corrections in this frame
and modifying it along the way. As the guess loop gradually approaches the
periodic orbit by repeatedly solving Eq. (8)) , modifications of the
coordinate system are made with the scheme illustrated in Fig. 2. In each
modification, each local coordinate frame moves along the loop for $n_{1}$
steps (roughly $5\sim 15$ steps in the current work). If the total number of
modifications during the whole search is $n_{2}$, in order to largely cover
the expanding directions, $n_{1}n_{2}\Delta t\geq 1/\lambda_{L}$ should be
satisfied, where $\lambda_{L}$ is some kind of average Lyapunov exponent. In
the following, we find that $n_{1}\times n_{2}>N/2$ could do the job quite
well, where $N$ is the number of lattice points depicting the periodic orbit.

We start the cycle search by numerical integration of the dynamical system
under investigation and locate possible orbit segments that nearly repeat
themselves - the recurring segments. Numerical experiments reveal regions
where a trajectory spends most of its life, giving us the first hunch where to
initialize a loop. After an extensive search for the recurring orbits, to
start the variational scheme, we first take a spatial FFT of one such segment
and keep only the lowest wavenumber components, of which an inverse Fourier
transform back to the phase space yields a smooth loop , which will be used as
our initial guess. In the following, we will apply the reduced scheme of the
variational method discussed above to several examples.

##  IV Application and discussion

In this section, we will use some examples to show the validity of the reduced
variational method. The convergence criterion in the following is

$Err.=\frac{1}{N}\sqrt{\sum\limits_{i=1}^{N}(\tilde{v}_{i}-\lambda
v_{i})^{2}}\leq 10^{-6}.$ (15)

###  IV.1 Navier-Stokes equation with periodic boundary conditions

The incompressible Navier-stokes equation on the torus
$T^{2}=[0,2\pi]\times[0,2\pi]$ is written as[19, 20]

$\displaystyle\frac{\partial\bm{u}}{\partial t}+(\bm{u}\cdot\nabla)\bm{u}=$
$\displaystyle-\nabla p+f+\nu\Delta\bm{u},$ (16) $\displaystyle div~{}\bm{u}=$
$\displaystyle 0,$ (17) $\displaystyle\int_{T^{2}}\bm{u}~{}dx=$ $\displaystyle
0,$ (18)

where $\bm{u}$ and $p$ is the velocity field and the pressure respectively,
and $f$ is the external periodic driving. Under certain conditions, Eq. (16)
could be expanded with $9$ Fourier modes[21], which reads

$\displaystyle\dot{x_{1}}=$
$\displaystyle-2x_{1}+4x_{3}x_{5}+4x_{2}x_{4}-4x_{6}x_{9}-4x_{7}x_{8},$
$\displaystyle\dot{x_{2}}=$ $\displaystyle-9x_{2}+3x_{1}x_{4},$
$\displaystyle\dot{x_{3}}=$ $\displaystyle-5x_{3}-x_{1}x_{5},$
$\displaystyle\dot{x_{4}}=$
$\displaystyle-5x_{4}-7x_{1}x_{2}+\frac{9\sqrt{5}}{5}x_{1}x_{7}+Re,$
$\displaystyle\dot{x_{5}}=$ $\displaystyle-
x_{5}-\sqrt{5}x_{1}x_{6}-3x_{1}x_{3},$ (19) $\displaystyle\dot{x_{6}}=$
$\displaystyle-x_{6}+3x_{1}x_{9}+\sqrt{5}x_{1}x_{5},$
$\displaystyle\dot{x_{7}}=$
$\displaystyle-5x_{7}-\frac{9\sqrt{5}}{5}x_{1}x_{4}+7x_{1}x_{8},$
$\displaystyle\dot{x_{8}}=$ $\displaystyle-9x_{8}-3x_{1}x_{7},$
$\displaystyle\dot{x_{9}}=$ $\displaystyle-5x_{9}+x_{1}x_{6},$

where $R_{e}=58$ is the Reynolds number.

Figure 3: The initial guess and the detected periodic orbit of a reduced version Eq. (IV.1) of the Navier-Stokes equation. (a) and (b) the projections of the initial guess loop in the $x_{1}-x_{4}$ and the $x_{7}-x_{8}$ plane, respectively. (c) and (d) for the periodic orbit obtained by the reduced variational method, with the parameter $\lambda=0.2216$ (indicating the period is 1.3924), $R_{e}=58$, the number of grid points $N=600$ and the dimension of the reduced coordinate frame $n=4$. Table 1: The Lyapunov exponents of the dynamical system described by Eq. (IV.1) for $Re=58$.  $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$ | $\lambda_{6}$ | $\lambda_{7}$ | $\lambda_{8}$ | $\lambda_{9}$
---|---|---|---|---|---|---|---|---
0.61292 | 0.00143 | -0.00069 | -1.81815 | -3.74089 | -6.52364 | -7.58614 | -8.61472 | -14.33010

For the system described by Eq. (IV.1), we calculated numerically its Lyapunov
exponents which are displayed in Table 1. It can be seen that the first two
values are greater than $0$. According to the previous discussion in
Sec.III.1, the dimension of the local coordinate system is set to $n=4$ (the
value of $\lambda_{3}$ and $\lambda_{4}$ are closest to $0$, which correspond
to the orbit direction and the translation symmetry of Eq. (IV.1)). The
results by the reduced method are displayed in Fig. 3 in which (a) and (b)
depict the projections of the initial guess loop constructed by the FFT, while
(c) and (d) plot the projection of the final cycle. The initial guess loop
does not seem to have any symmetry, but the cycle is obviously symmetric. As
can be seen from diagrams (a) and (c), the initial loop roughly possesses the
topology of the periodic orbit it is about to approach. Therefore, even if the
initial conditions are not very good, our method still drives the guess loop
to a periodic orbit.

Figure 4:  The convergence of the reduced variational method for Eq. (IV.1)
with $N=1000$, $900$, $800$, $700$, $600$, $Re=58$ and the dimension of
reduced space is $n=4$

Here we further investigate the influence of the number of lattice points on
the convergence. The same periodic orbit (as shown in Fig. 3) is used as the
target with different numbers of lattice points
$N=1000,~{}900,~{}800,~{}700,~{}600$ and the accuracy of the initial loop
being kept the same. With the guess loop iterated for $100$ times, the
convergence of our reduced method is shown in Fig. 4. It is clear that in the
first $40$ iterations, the convergence are similar among different
representations, but becomes different in the ensuing iterations. The larger
the number of orbit points, the better the convergence. The fast convergence
at the first stage stems from the initial exponential decay of the error in
the projected subspace and thus is similar for different discretization.
However, the relatively slow convergence at the second stage is enabled by the
modifications of the coordinate frames which proceed distinctively for
different $N$’s.

###  IV.2 The Coupled Lorenz Systems

In the flowing we test our reduced method in case of the coupled Lorenz
systems[22] $\sum_{i}:\\{\dot{x}_{i}(t),\dot{y}_{i}(t),\dot{z}_{i}(t)\\}$,
described by the following equations

$\displaystyle\dot{x_{i}}=$ $\displaystyle\sigma(y_{i}-x_{i}),$
$\displaystyle\dot{y_{i}}=$ $\displaystyle x_{i}(\rho-
z_{i})-y_{i}+k_{i}y_{i+1},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(i=1,2,3,...,\bar{n})$
(20) $\displaystyle\dot{z_{i}}=$ $\displaystyle x_{i}y_{i}-\mu z_{i},$

where the parameters $\sigma=10$, $\rho=28$, $\mu=8/3$ with a periodic
boundary condition: $y_{\bar{n}+1}=y_{1}$. Here, $\bar{n}$ is the number of
the coupled lorenz systems and $k_{i}$ is the coupling strength. In the
uncoupled case $k_{i}=0$, all the subsystem perform chaotic motions on the
well-known Lorenz attractor. Linear coupling is imposed via the $y_{i+1}$
components of the adjacent subsystem, controlled by $k_{i}$. Here, we take the
case of $\bar{n}=33$ as an example, where the total dimension of the whole
system is $d=99$. In order to make the coupling of the system non-uniform, we
set the coupling strength as $k_{i}=[0.5+0.5\frac{(i-1)}{(\bar{n}-1)})]k_{0}$
(with $k_{0}=0\sim 1$).

Table 2: The Lyapunov exponents of the $33$ coupled Lorentz systems with $k_{0}=1$. Here we only show the first 35 Lyapunov exponents. $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$ | $\lambda_{6}$ | $\lambda_{7}$ | $\lambda_{8}$ | $\lambda_{9}$ | $\lambda_{10}$ | $\lambda_{11}$ | $\lambda_{12}$
---|---|---|---|---|---|---|---|---|---|---|---
1.3726 | 1.2824 | 1.2119 | 1.1894 | 1.1337 | 1.1154 | 1.0835 | 1.0109 | 0.9931 | 0.9299 | 0.8913 | 0.8822
$\lambda_{13}$ | $\lambda_{14}$ | $\lambda_{15}$ | $\lambda_{16}$ | $\lambda_{17}$ | $\lambda_{18}$ | $\lambda_{19}$ | $\lambda_{20}$ | $\lambda_{21}$ | $\lambda_{22}$ | $\lambda_{23}$ | $\lambda_{24}$
0.8017 | 0.8005 | 0.7710 | 0.6820 | 0.6550 | 0.6428 | 0.5747 | 0.5522 | 0.5258 | 0.4586 | 0.4470 | 0.3925
$\lambda_{25}$ | $\lambda_{26}$ | $\lambda_{27}$ | $\lambda_{28}$ | $\lambda_{29}$ | $\lambda_{30}$ | $\lambda_{31}$ | $\lambda_{32}$ | $\lambda_{33}$ | $\lambda_{34}$ | $\lambda_{35}$ | …
0.3609 | 0.3086 | 0.2899 | 0.2871 | 0.2331 | 0.2098 | 0.1945 | 0.1635 | 0.1161 | 0.0176 | 0.01184 | …

Figure 5:  Periodic orbits of the coupled lorentz systems. (a) and (b), the
periodic orbits for the first and the 33 _rd_ Lorenz system with the coupling
$k_{i}=0$ and $T=1.5997$; (c) the linear relation of the variable $y_{1}$ and
$y_{33}$ during the evolution; (d) and (e), the periodic orbit obtained by the
reduced method with the coupling $k_{0}=1$, $T=1.5915$, with $N=600$ and
$n=34$; (f) the nonlinear dependence of $y_{1}$ on $y_{33}$.

For such a system, we may employ an alternative way of constructing an
interesting initial loop. Since the system consists of $33$ coupled lorentz
oscillators, we may use periodic orbits of the original lorentz system to
construct the initial loop. For example, we take the same cycle for each
lorenz system (as shown in Fig.5 (a)$\sim$(c), the periodic orbit is
symmetrical with respect to $x=0$). When the coupling strength is $k_{i}=0$,
the initial guess is already a true periodic orbit of the whole system.
Recursively, we search a nearby new cycle upon gradually increasing the
coupling constant, starting with the already determined cycle at the previous
coupling. This provides an alternative way to start the search and a family of
cycles are obtained when $k_{0}$ goes from $0$ to $1$. We have numerically
calculated the Lyapunov exponents displayed in the Table 2, and found that the
values of the first $33$ Lyapunov exponents are positive and away from zero
corresponding to the main stretching directions. Such a result is also easy to
understand, because there is an expanding direction in each Lorentz system,
therefore the $33$ coupled Lorentz systems should have $33$ unstable
directions as long as the coupling is small. According to the previous
discussion in Sec.III.1, the dimension of the reduced space is set as $n=34$.

The final result is shown in Figs.5 (d) $\sim$ (f). By comparing the graphs in
(d) and (e), we see that the periodic orbit of each individual Lorentz system
does not seem to change much in the $x-y$ plane. However, the relation between
corresponding variables in different subsystems exhibit nontrivial features
(see Figs. 5 (f)). In this example, from the construction of the initial loop
to the search of periodic orbit, the whole process is completed in local
subspaces without evolution in the full phase space, which certainly cuts down
the computation complexity and provides a good example to deal with high-
dimensional systems with similar spatial structures. As shown in Fig.6, the
algorithm converges faster initially and slows down when the guess loop is
about to approach the true periodic orbit for the same reason given in the
previous example. But on the whole, the convergence is good and remains
exponential.

In previous example, the unstable directions are quite few so that the
dimension of the reduced local coordinate subspace is also low, the
convergence may be expected. Nevertheless, in the current example, there are
many unstable directions but the convergence is still good with a proper
construction of the local subspaces. As a matter of fact, the local dimension
can be further reduced even down to $n=9$ and the method is still effective
though with a slow convergence. How is it possible? During the evolution, all
the local coordinate frames are essentially different. Each modification
rotates these frames to new dimensions and the variational method will bring
down errors in these directions. Finally, when all the directions are covered,
a cycle may be identified. Nevertheless, the number of the included dimensions
could not be too small to make the algorithm unstable. We still have no idea
of what the necessary number of dimension is.

Figure 6:  The convergence of the reduced variational method for the coupled
lorenz system with $k_{0}=1$ and the dimension of the reduced space is $n=34$

###  IV.3 The Kuramoto-Sivashinsky system

Figure 7:  The Kuramoto-Sivashinsky system in a spatiotemporally chaotic
regime (viscosity parameter $\nu=0.015$, $d=16$ Fourier mode truncation). (a)
and (b), the projection of the initial guess loop in the $x_{1}-x_{2}$ plane
and the $x_{14}-x_{16}$ plane, respectively. (c) and (d), the projection of
the resulting periodic orbit with $\lambda=0.0938$, $T=0.5892$, $N=856$ and
$n=5$.

Here, we directly apply the current scheme to a nonlinear partial differential
equation - the Kuramoto-Sivashinsky equation (KSe) to check the influcence of
high dimensions and possible complications arising from possible stiffness due
to the existence of inertial manifold. The KSe was independently proposed by
Kuramoto and Tsuzuki in the study of phase turbulence in reaction-diffusion
systems [23] and Sivashinsky in the study of flame combustion propagation
models [24]. The KSe is an intrinsic phase equation describing the slow change
of vibrations in a spatially extended system and also is an ideal model for
studying high-dimensional systems. The KSe in one-dimensional space reads

$u_{t}=(u^{2})_{x}-u_{xx}-\nu u_{xxxx},$ (21)

where the subscript $x$ and $t$ represent the partial derivative of $u$ with
respect to $x$ or $t$, with $t\geq 0$ and $x\in[0,2\pi]$. The $\nu$ is a
viscous damping parameter which is the dissipation rate of the system, and the
first term on the right hand side of Eq. (21) is a nonlinear convection term,
which induces interaction between different spatial modes and energy transfer
from low wave-number modes to high ones. As $\nu$ decrease, the system
undergoes a series of bifurcations, leading to increasingly turbulent
dynamics.

If we choose to study only the odd solutions $u(t,-x)=-u(t,x)$ with the
periodic boundary condition $u(t,x+2\pi)=u(t,x)$, the state variable $u(x,t)$
can be expanded in spatial Fourier series [25],

$u(t,x)=i\sum\limits_{k=-\infty}^{\infty}a_{k}(t)e^{ikx},$ (22)

where $a_{-k}=-a_{k}\in\mathbb{R}$. In terms of the Fourier components, Eq.
(21) is rewritten as an infinite ladder of ODEs:

$\dot{a}_{k}=(k^{2}-\nu
k^{4})a_{k}-k\sum\limits_{k=-\infty}^{\infty}a_{m}a_{k-m}.$ (23)

In numerical computation, a Galerkin truncation of the Fourier series has to
be applied. As the amplitude of $a_{k}$ decreases with the increase of $k$
near the strange attractor, the high wavenumber modes in the asymptotic regime
are negligible. Thus we can reduce the equation to a finite but large number
of ODEs by the Galerkin truncation. In the current exploration, we work with
two different truncations $d=16$ and $64$ in the following. However, in Eq.
(23), the fourth power of $k$ arising from the dissipation could be quite
large for large $k$, which makes the system stiff and brings trouble to
numerical computation.

Table 3: The Lyapunov exponents of The Kuramoto-Sivashinsky system with $d=16,~{}\nu=0.015$ and $d=64,~{}\nu=0.0049$. Here we only show the first 8 largest Lyapunov exponents.  $d=16$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$ | $\lambda_{6}$ | $\lambda_{7}$ | $\lambda_{8}$
---|---|---|---|---|---|---|---|---
12.8759 | 3.9500 | -0.0092 | -2.8395 | -8.1252 | -13.6478 | -18.4199 | -23.5989
$d=64$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$ | $\lambda_{6}$ | $\lambda_{7}$ | $\lambda_{8}$
12.9183 | 7.8012 | 3.7726 | 0.3122 | -0.0013 | -3.2319 | -8.1940 | -20.0568

The Lyapunov exponents of the system are computed and filled in Table. 3. We
show the values of the first eight Lyapunov exponents at $d=16$ ($\nu=0.015$)
and $64$ ($\nu=0.0049$), where the numbers of positive values are $2$ and $4$.
The exponents $\lambda_{3}=-0.0092$ and $\lambda_{5}=-0.0013$ correspond to
the orbit direction of the system with $d=16$ and $d=64$, respectively.
According to the previous discussion in Sec.III.1, we select $n=5$ and $n=8$
directions to form the local coordinate frames along the loop, including the
velocity direction.

Figure 8:  The Kuramoto-Sivashinsky system in a spatiotemporally turbulent
regime (viscosity parameter $\nu=0.0049$, $d=64$ Fourier mode truncation). (a)
and (b), the projection of the initial guess loop on the $x_{11}-x_{13}$ plane
and the $x_{56}-x_{63}$ plane, respectively. (c) and (d), the projection of
resulting periodic orbit with $\lambda=0.0263$, $T=0.1649$, $N=467$ (the
number of grid points) and $n=8$ (The dimension of the reduced coordinate
frame).

As shown in Figs. 7 and 8, the projections of the determined periodic orbits
and the corresponding initial guess loops of with $d=16$ and $64$ are shown.
It is found that the amplitudes of high-wavenumber modes are much smaller than
the low ones ($\bm{a}_{11}/\bm{a}_{63}\sim 10^{4}$).

Nevertheless, the high-wavenumber modes have intricate features, which makes
them sensitive to the initial guess, especially when it is rough. This is not
surprising considering the stiffness problem mentioned after Eq. (23).
Therefore, extra care should be exercised during the modification of high-
wavenumber modes the corrections of which could be multiplied by someweight
ranging from $1$ to $0.1$ to slow down the variation. In addition, it is not
difficult to find that the topological structures of the high-wavenumber modes
change a great deal during loop evolution when the guess loop is not good,
which prevents a smooth convergence and possibly induces numerical
instability. With further increase of the dimensions of the system, this
problem becomes prominent (such as in the case of $d=64$).

Figure 9:  The convergence of the reduced variational method for the Kuramoto-
Sivashinsky system with $\nu_{0}=0.0049,~{}d=64$ Fourier mode truncation and
the dimension of reduced space is $8$.

In order to promote the robustness of the reduced method the average position
of the neighboring points on the loop is used as the new value at $s_{i}$
(namely, $\bm{X}(s_{i})\longrightarrow[\bm{X}(s_{i})+\bm{X}(s_{i+1})]/2$)
after each correction. Despite all the details dealing with a bad initial
guess, in the Kuramoto-Sivashinsky system with $d=64$ Fourier mode truncation,
for a subspace with dimension $n=8$, our reduced method has good convergence.
As shown in Fig. 9, when the initial guess loop is close to a periodic orbit
(entering the linearization neighborhood of the orbit), the error function
Err. will decreases exponentially.

Figure 10:  Level plot of the space-time evolution u(x,t) in the (t, x) plan.
(a) and (c), the initial guesses in Fig. 7 and Fig. 8 respectively. (b) and
(d), the periodic orbits depicted in Fig. 7 (with the cycle period $T=0.5892$
) and Fig. 8 (a repeat followed after the cycle period $T=0.1649$ ).

The space-time evolutions of $u(x,t)$, the unstable spatiotemporally solutions
corresponding to the orbits in Fig. 7 and Fig. 8, are plotted in Fig. 10. As
$u(x,t)$ is antisymmetric in $[-\pi,\pi]$, it suffices to display the
solutions in the $[0,\pi]$ interval. There is little difference between the
spatiotemporal profiles in Fig. 10 (a) and (b), which indicates that the guess
loop is close to the periodic orbit. While big difference between (c) and (d)
indicates that the initial guess is bad. A comparison of (b) and (d) reveals
that the structure displayed in (d) has more features than that in (b) which
indicates that the higher the dimension of the attractor, the finer the space-
time evolution of the $u(x,t)$.

Figure 11:  The estimated storage and time cost required by different
variational approaches. (a) the storage resource (SR) required by the original
variational method (VM, black-quare curve) and by the reduced method (RM,
black-dot curve). The number of lattice points is $10^{3}$. (b) the estimated
time required by the RM (black-dot curve) and the VM (black-quare curve). The
red-circle marks the ratio ($t_{VM}/t_{RM}$). The number of grid points is
$746,~{}856,~{}856,$ and $467$ for systems with dimensions $d=9,~{}16,~{}32$
and $64$. It should be noted here that the system corresponding to $d=9$ is
the previous nine-mode Lorentz system.

Often than not, the original variational method (VM) consumes a lot of storage
resource (SR) and computer time which will become more and more striking with
the increase of system dimension, but the reduced variational method (RM)
introduced in the current manuscript partly solves the problem. As shown in
Fig. 11, we compared the required storage and computing time for the VM and
the RM. In Fig. 11 (a), the black-quare curve depicts the storage resource
required by the VM, which increases exponentially with the increase of system
dimension $d$, while the black-dot curve for the RM grows only linearly. The
difference between them becomes more and more pronounced with the increase of
system dimension $d$ (the ratio of storage resources required is
$SR_{VM}/SR_{RM}\sim 6$ for $d=64$). In Fig. 11 (b), the black-square curve is
the time required by the VM, and the black-dot curve by the RM. Because for
different system dimensions, the number of lattice points is different, it is
difficult to compare directly. A quantity $t_{VM}/t_{RM}$ (that is, the ratio
of the time $t_{VM}$ required by VM to the time $t_{RM}$ required by RM)is
defined and plotted with the red-circle curve in Fig. 11 (b). It is clear that
with the increase of system dimension, the ratio $t_{VM}/t_{RM}$ becomes
larger and larger, and reaches $t_{VM}/t_{RM}\sim 28$ at $d=64$. Certainly,
when the dimension of the system is not high, the advantage of the RM is not
so prominent. Henceforth, the caveat that lies in the variational method is
essentially filled in the reduced scheme and we expect its fruitful
application to high dimensional systems.

##  V SUMMARY

The variational method proposed in Ref. [13] is very robust numerically and
hence widely used when the system dimension is not high. However, in
application to complex systems with high dimensions, its bottleneck of
computation complexity is becoming prominent. That is, in solving the
associated matrix equation, the demand for storage and computing time
increases sharply with the increase of the phase space dimension. In order to
break the bottleneck, we proposed a reduction scheme in the current work. The
main idea is to reduce the high-dimensional coordinate frame at each loop
point to a low-dimensional one based on the existence of inertial manifold
widely observed in various spatially extended systems [26]. When the guess
loop resides in the neighborhood of a periodic orbit, to enable the
convergence, the unstable directions need to be carefully controlled, while
the stable directions drag the loop to the periodic orbit during forward
evolution. Based this consideration, we give a reduced version of the
variational method to take care of all the unstable directions but implement
the natural evolution in the stable ones. In addition, a pseudo-evolution
equation along the loop and its improved version are proposed to evolve the
local coordinate systems. An alternating evolution of the loop and the local
coordinate frames is practiced and observed to lead to robust convergence.

In all examples, a rough criterion for determining the dimension $n$ of the
reduced space is verified and should be effective for general systems that $n$
is greater than the number of positive Lyapunov exponents of the system. The
velocity direction needs to be included in the construction of the local
coordinate system, because in the early iteration of the guess loop, together
with the expanding directions it enables a fast convergence. However, other
directions may take the lead in subsequent iterations, which is also the
reason why there is an inflection point on the convergence curve. Therefore,
in order to make the algorithm stable and converging well, the alternation
between the loop evolution ($10\sim 20$ iterations) and the modification of
local coordinates (about $5\sim 15$ steps) is carried out and found to be
effective. Away from the periodic orbit, the reduced scheme could become
unstable just like the original variational method. But through an averaging
process, the stability is restored quickly as discussed in Section IV.3 in its
application to the $64$-dimensional Kuramoto-Sivashinsky system. When the
guess loop enters the linearization neighborhood of a periodic orbits, the
algorithm shows very good exponential convergence.

In the new reduced scheme, the evolution of local coordinate systems is based
on pseudo-evolution. If the initial loop is poor, it is hard to pin down the
stretching directions so that the convergence is quite slow, which makes the
algorithm inefficient. Better schemes should be designed to utilize and
benefit from possible low-dimensional structures. Nevertheless, the reduced
scheme shows superior stability and efficiency for searching periodic orbits
in high-dimensional systems, and seems to provide a new tool to explore
complex system dynamics dominated by recurring patterns.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
under Grants No. 11775035, and also by the Fundamental Research Funds for the
Central Universities with Contract No.2019XD-A10.

## References

  * Kawahara and Kida [2001a] G. Kawahara and S. Kida, “Periodic motion embedded in plane couette turbulence: regeneration cycle and burst,” J. Fluid Mech. 449, 291–300 (2001a).
  * Kato and Yamada [2003] S. Kato and M. Yamada, “Unstable periodic solutions embedded in a shell model turbulence,” Phys. Rev. E 68, 025302 (2003).
  * Gutzwiller [1990] M. C. Gutzwiller,  _Chaos in Classical and Quantum Mechanics_ (Springer, New York, 1990).
  * Cvitanović [1988] P. Cvitanović, “Invariant measurement of strange sets in terms of cycles,” Phys. Rev. Lett. 61, 2729–2732 (1988).
  * Artuso, Aurell, and Cvitanovic [1990] R. Artuso, E. Aurell, and P. Cvitanovic, “Recycling of strange sets: I. cycle expansions,” Nonlinearity 3, 325–359 (1990).
  * Lan [2010] Y. Lan, “Cycle expansions: From maps to turbulence,” Commun. Nonl. Sci. 15, 502–526 (2010).
  * Auerbach  _et al._ [1987] D. Auerbach, P. Cvitanović, J.-P. Eckmann, G. Gunaratne, and I. Procaccia, “Exploring chaotic motion through periodic orbits,” Phys. Rev. Lett. 58, 2387–2389 (1987).
  * Mestel and Percival [1987] B. Mestel and I. Percival, “Newton method for highly unstable orbits,” Physica D 24, 172–178 (1987).
  * Baranger, Davies, and Mahoney [1988] M. Baranger, K. Davies, and J. Mahoney, “The calculation of periodic trajectories,” Ann. Phys. 186, 95–110 (1988).
  * Parker and Chua [1989] T. S. Parker and L. Chua,  _Practical Numerical Algorithms for Chaotic Systems_ (Springer, New York, 1989).
  * Farantos [1995] S. C. Farantos, “Methods for locating periodic orbits in highly unstable systems,” J. Mol. Struct. Theochem 341, 91–100 (1995).
  * Deane and Marsh [2006] J. Deane and L. Marsh, “Unstable periodic orbit detection for odes with periodic forcing,” Phys. Lett. A 359, 555–558 (2006).
  * Lan and Cvitanovic [2004] Y. Lan and P. Cvitanovic, “Variational method for finding periodic orbits in a general flow,” Phys. Rev. E 69, 016217 (2004).
  * Saiki [2007] Y. Saiki, “Numerical detection of unstable periodic orbits in continuous-time dynamical systems with chaotic behaviors,” Nonlinear Proc. Geoph. 14, 615–620 (2007).
  * Kawahara and Kida [2001b] G. Kawahara and S. Kida, “Periodic motion embedded in plane couette turbulence: regeneration cycle and burst,” J. Fluid Mech. 449, 291–300 (2001b).
  * Zhou, Ren, and Weinan [2008] X. Zhou, W. Ren, and E. Weinan, “Adaptive minimum action method for the study of rare events,” J. Chem. Phys. 128, 104111 (2008).
  * Dong and Lan [2014] C. Dong and Y. Lan, “A variational approach to connecting orbits in nonlinear dynamical systems,” Phys. Lett. A 378, 705–712 (2014).
  * Wang, Wang, and Lan [2018] D. Wang, P. Wang, and Y. Lan, “Accelerated variational approach for searching cycles,” Phys. Rev. E 98, 042204 (2018).
  * Boldrighini and Franceschini [1979] C. Boldrighini and V. Franceschini, “A five-dimensional truncation of the plane incompressible navier-stokes equations,” Commun. math. Phys. 64, 159–170 (1979).
  * Franceschini and Tebaldi [1981] V. Franceschini and C. Tebaldi, “A seven-mode truncation of the plane incompressible navier-stokes equations,” J. Stat. Phys. 25, 397–417 (1981).
  * Yan [2014] G. Yan, “The dynamical behavior and the numerical simulation of a nine-mode lorenz equations,” Acta Math. Appl. Sin. 37, 507–515 (2014).
  * Frenzel and Pompe [2007] S. Frenzel and B. Pompe, “Partial mutual information for coupling analysis of multivariate time series,” Phys. Rev. Lett. 99, 204101 (2007).
  * Kuramoto and Tsuzuki [1976] Y. Kuramoto and T. Tsuzuki, “Persistent propagation of concentration waves in dissipative media far from thermal equilibrium,” Progr. Theor. Phys. 55, 365–369 (1976).
  * Michelson and Sivashinsky [1977] D. M. Michelson and G. I. Sivashinsky, “Nonlinear analysis of hydrodynamic instability in laminar flames—ii. numerical experiments,” Acta Astr. 4, 1207–1221 (1977).
  * Christiansen, Cvitanovi$\acute{\text{c}}$, and Putkaradze [1997] F. Christiansen, P. Cvitanovi$\acute{\text{c}}$, and V. Putkaradze, “Spatiotemporal chaos in terms of unstable recurrent patterns,” Nonlinearity 10, 55–70(16) (1997).
  * Temam [1988] R. Temam,  _Infinite-Dimensional Dynamical Systems in Mechanics and Physics_ (Springer-Verlag, New York, 1988).


AGC249525 | 99.26(25) | 2.09 (1.93-2.12) | 14:17:53.9 | +17:32:46.0 | 6.80${}^{+0.01}_{-0.03}$ | -5.99${}^{+0.07}_{-0.03}$ | -6.87${}^{+0.07}_{-0.03}$ | 0.96$\pm$0.03 | 0.92$\pm$0.02 | 4.77${}^{+0.34}_{-0.38}$ | 5.09${}^{+0.01}_{-0.05}$ | 106.40$\pm 2.15$ | 50.84$\pm 0.54$