# Balancing Cooperativeness and Adaptiveness in the (Noisy) Iterated
Prisoner’s Dilemma

Adrian Hutter
<EMAIL_ADDRESS>
###### Abstract

Ever since Axelrod’s seminal work, tournaments served as the main benchmark
for evaluating strategies in the Iterated Prisoner’s Dilemma (IPD). In this
work, we first introduce a strategy for the IPD which outperforms previous
tournament champions when evaluated against the 239 strategies in the Axelrod
library, at noise levels in the IPD ranging from 0% to 10%. The basic idea
behind our strategy is to start playing a version of tit-for-tat which
forgives unprovoked defections if their rate is not significantly above the
noise level, while building a (memory-1) model of the opponent; then switch to
a strategy which is optimally adapted to the model of the opponent. We then
argue that the above strategy (like other prominent strategies) lacks a couple
of desirable properties which are not well tested for by tournaments, but
which will be relevant in other contexts: we want our strategy to be _self-
cooperating_ , i.e., cooperate with a clone with high probability, even at
high noise levels; and we want it to be _cooperation-inducing_ , i.e., optimal
play against it should entail cooperating with high probability. We show that
we can guarantee these properties, at a modest cost in tournament performance,
by reverting from the strategy adapted to the opponent to the forgiving tit-
for-tat strategy under suitable conditions.

##  1 Introduction

Direct reciprocity is one of the main mechanisms explaining the emergence of
cooperation between self-interested individuals (Nowak,, 2006), and the
Iterated Prisoner’s Dilemma (IPD) is the primary model used for the study of
direct reciprocity. Famously, tit-for-tat (TFT) won both of Robert Axelrod’s
seminal IPD tournaments (Axelrod,, 1984). In any realistic situation
resembling an IPD, intended actions can be executed imperfectly or otherwise
have unintended consequences. The simplest way of modelling this in the IPD is
to assume that each action has the opposite of the intended effect with some
probability $p_{\text{noise}}$. It was soon recognized that the noisy IPD
introduces fundamentally new challenges. In particular, for TFT a single noise
event ends mutual cooperation (at least till the next noise event occurs),
leading to alternating rounds of C/D, D/C instead (Axelrod and Hamilton,,
1981). This motivated the development of alternatives to TFT which are able to
cooperate more robustly in the presence of noise, including tit-for-two-tats
(Axelrod,, 1984), Generous TFT (Molander,, 1985), Pavlov (aka. Win-Stay Lose-
Shift) (Nowak and Sigmund,, 1993), and Contrite TFT (Wu and Axelrod,, 1995).
Wu and Axelrod, (1995) evaluated the latter three strategies against all 63
strategies submitted to Axelrod’s second tournament (Axelrod,, 1984) (using
$p_{\text{noise}}=1\%$), finding Generous and Contrite TFT to perform
strongly, and Pavlov to perform poorly.

In celebration of the 20th anniversary of Axelrod’s seminal work, Kendall et
al., (2007) organized two IPD tournaments (both noise-free and noisy, with
$p_{\text{noise}}=10\%$) in 2004 and 2005. We briefly describe two of the
best-performing strategies, both of which are designed specifically to account
for noise. Omega TFT (Slany et al.,, 2007) modifies TFT in two ways: if it
detects a deadlock loop (alternating rounds of C/D and D/C, which are
characteristic for TFT) it attempts to break the loop by playing C twice; and
it reverts to playing D for the remainder of the game if a measure of
randomness of the opponent’s play exceeds a threshold. DBS (Au and Nau,, 2007)
models the opponent as a memory-1 strategy.111I.e., the opponent is described
by the probabilities of playing C after the 4 possible states
$\\{\text{CC},\text{CD},\text{DC},\text{DD}\\}$. It attempts to describe the
opponent using deterministic rules, and ignores occasional violations of these
rules (which might be due to noise). It then optimizes its move against the
model of the opponent using tree search to depth 5.

More recently, Harper et al., (2017) trained various model types (lookup
tables, artificial neural networks, finite state machines, hidden Markov
models), in both noisy and noise-free environments, against a large zoo of
strategies from previous literature. They then ran a noisy
($p_{\text{noise}}=5\%$) and a noise-free tournament of 176 strategies,
including the trained strategies and all strategies mentioned so far.
Remarkably, DBS was the best-performing human-designed strategy in both
tournaments, ranking 1st in the noisy tournament and 12th in the noise-free
tournament. Omega TFT ranked not far behind DBS (8th in the noisy tournament
and 15th in the noisy one). The first ranked strategy in the noise-free
tournament was EvolvedLookerUp2_2_2, a deterministic lookup table, which bases
its decision on the first 2 actions of the opponent and the past 2 actions of
itself and the opponent. All strategies participating in these tournaments
(and more) are available in the Axelrod library (Knight et al.,, 2016), which
we will use for evaluation purposes in the present work.

The goal of the present work is two-fold. Firstly, we want to create a
strategy which is both able to robustly cooperate (like the forgiving variants
of TFT) and highly adaptive (by building a simple model of the opponent and
responding optimally to that, like DBS). These goals cannot be pursued at the
same time, so we need rules for switching between the cooperative and the
adaptive state. Sec. 2 introduces such a strategy and evaluates it against all
strategies in the Axelrod library. Secondly, we discuss properties which are
desirable for a strategy to have, beyond doing well in tournaments. In Sec. 3,
we argue that a strategy should ideally be _self-cooperating_ (cooperate with
a clone) and _cooperation-inducing_ (incentivize the opponent to cooperate).
We then show in Sec. 4 how to make the previously introduced strategy self-
cooperating and cooperation-inducing by introducing additional rules for
switching between the cooperative and the adaptive state.

##  2 Cooperate, then adapt to the opponent

The strategy introduced in this section, which we call _CooperateISO_ ,
comprises two sub-strategies: a highly cooperative one, and a purely adaptive
one, together with a rule for switching from the former to the latter. The
cooperative strategy, _Longterm TFT_ focuses on allowing robust cooperation in
the presence of noise without being exploitable. The adaptive strategy, which
we call _ISO_ (Infinite Sum Optimizer), builds a memory-1 model of the
opponent and responds optimally to that. The basic idea is illustrated in Fig.
1, and motivated and explained in more detail in the following subsections.

Figure 1: The strategy introduced in this work comprises two sub-strategies,
Longterm TFT (focused on allowing robust cooperation in the presence of noise)
and ISO (focused on playing optimally against simple opponents), together with
rules for switching between these. Sec. 2 introduces the two sub-strategies
and the rule for switching from Longterm TFT to ISO; Sec. 4 introduces two
rules for switching back, which are motivated in Sec. 3.

In the following, we use the standard notation for referring to the possible
payoffs in the PD, $T>R>P>S$, and use the conventional values $T=5$, $R=3$,
$P=1$, and $S=0$.

###  2.1 Longterm TFT: a maximally forgiving version of TFT

Multiple existing strategies feature mechanisms for dealing with the effects
of noise (e.g. Molander, (1985); Nowak and Sigmund, (1993); Wu and Axelrod,
(1995); Slany et al., (2007); Au and Nau, (2007)). In this subsection, we
introduce a strategy which, like Generous TFT (Molander,, 1985), is designed
to be maximally forgiving (of defections which can be due to noise) while
still incentivizing the opponent to cooperate. However, while Generous TFT is
a memory-1 strategy, our strategy takes (aggregate statistics of) the entire
history of play into account; we thus call this strategy _Longterm TFT_.

Longterm TFT starts by playing standard TFT. It switches to always cooperating
if the history of play is compatible (given the noise model) with the opponent
always rewarding cooperation with cooperation. More precisely, let $N_{C}$ be
the number of times our agent has cooperated so far, not taking our agent’s
most recenct action into account (since the opponent has not yet reacted to
that). (Note that we care about the _actual_ action of our agent here, taking
the effects of noise into account, not the _intended_ action.) Let
$N_{C\rightarrow D}$ be the number of times the opponent defected after our
agent cooperated. Define the statistic

$\displaystyle z=\frac{N_{C\rightarrow
D}-p_{\text{noise}}N_{C}}{\max\\{1,\sqrt{p_{\text{noise}}(1-p_{\text{noise}})N_{C}}\\}}\
.$ (1)

Under the null hypothesis that the opponent always rewards cooperation with
cooperation, $N_{C\rightarrow D}$ is binomially distributed with mean
$p_{\text{noise}}N_{C}$ and variance
$p_{\text{noise}}(1-p_{\text{noise}})N_{C}$. As a consequence of the Central
Limit Theorem, $z$ will thus follow a standard normal distribution under the
null hypothesis, if $p_{\text{noise}}(1-p_{\text{noise}})N_{C}\geq 1$ and
$N_{C}$ is “large enough”. The statistic $z$ thus quantifies how confidently
we can reject the null hypothesis that the opponent always rewards cooperation
with cooperation.

Longterm TFT always cooperates if $N_{C}\geq 5$ and $z<2$.222No attempt was
made to optimize the precise threshold values, and similar thresholds in the
rest of this work. For $p_{\text{noise}}N_{C}<1$, i.e., if we expect less than
one unprovoked defection under the null hypothesis, the denominator in Eq. (1)
becomes $1$, and so Longterm TFT forgives at most 2 unprovoked defections.

During evaluation, we will assume that $p_{\text{noise}}$ is known at the
beginning of the IPD. If it is not known, it has to be estimated from
agreement between intended and actual actions. Since Longterm TFT always plays
identically to TFT while $N_{C}<5$, it does not need an estimate of
$p_{\text{noise}}$ during the first few rounds of the IPD.

By forgiving defections which are not significantly more frequent than one
would expect based on the noise model, Longterm TFT prevents unnecessary
chains of retaliation when playing against a cooperative but provocable
strategy like TFT and its variations. It is thus able to robustly cooperate
with such strategies (including itself) even in the presence of strong noise.

By being so forgiving, Longterm TFT invites occasional deliberate defections
from the opponent. However, the number of defections which are forgiven only
grows like $O(\sqrt{N_{C}})$, and so the _rate_ of defections which are
forgiven goes to zero like $O(1/\sqrt{N_{C}})$. The “free” defections which
Longterm TFT forgives are thus irrelevant in the steady state. In the steady
state, Longterm TFT always cooperates against all strategies which are at
least as cooperative as TFT: they may retaliate against Longterm TFT’s
(accidental) defections, but have to reward cooperation with cooperation.

###  2.2 ISO: playing optimally against memory-1 opponents

Many of the most prominent strategies in the IPD are memory-1, i.e., their
probability of cooperating on the next move only depends on the actions played
by both players in the previous round. Such strategies are thus described by
four probabilities.333 Plus a fifth probability describing the probability of
cooperating in the very first round of the IPD, which is irrelevant for our
purposes. DBS (Au and Nau,, 2007) builds a memory-1 model of the opponent and
then optimizes its own actions against that model. The strategy described in
this subsection, _ISO_ (Infinite Sum Optimizer), takes inspiration from DBS’
remarkably strong performance in previous tournaments. ISO is both a
simplification and a refinement of DBS.

DBS follows a complex algorithm for building its opponent model, trying to
establish deterministic rules (such as “the opponent always cooperates after
mutual cooperation”) while taking possible noise events into account. By
contrast, ISO’s opponent model is simply given by the (discounted) average
rate of cooperation of the opponent in the four different states. We use a
discount factor $\gamma_{\text{past}}=0.99$ in order to react more quickly to
changes in the opponent’s strategy than if a simple average were taken. In
order to make the averages well-defined from the beginning, we assume that we
have seen the opponent play the action dictated by TFT once for each of the
four states. We also clamp all empirical rates of cooperation to the interval
$[p_{\text{noise}},1-p_{\text{noise}}]$, since values outside of that interval
are not feasible under the noise model. Let $\vec{p}_{\text{opp}}^{\,(n)}$
denote our model of the opponent obtained this way, i.e., the 4-vector of
discounted average rates of cooperation in the four different states
$(CC,CD,DC,DD)$. We use the superscript $(n)$ to denote that these rates of
cooperation already include the effects of noise, and so might be different
from the _intended_ rates of cooperation.

While DBS uses tree search to depth 5 in order to optimize its action against
the opponent model, ISO optimizes the exact discounted sum of future payoffs.
Let $\vec{p}_{\text{self}}$ denote the 4-vector of probabilities of
cooperation of our agent, which we hope to optimize against
$\vec{p}_{\text{opp}}^{\,(n)}$. In order to calculate how well a candidate
$\vec{p}_{\text{self}}$ performs against a given
$\vec{p}_{\text{opp}}^{\,(n)}$, we first calculate
$\vec{p}_{\text{self}}^{\,(n)}$ by taking the effects of noise into account,
i.e., applying $p\mapsto(1-p_{\text{noise}})p+p_{\text{noise}}(1-p)$ to each
entry in $\vec{p}_{\text{self}}$. Iterated play between two memory-1
strategies (ISO and the model of the opponent) induces a Markov process of
order 1. We can then calculate the 4x4 transition matrix of this Markov
process,

$\displaystyle
T=\left(\vec{p}_{\text{self}}^{\,(n)}\odot\vec{p}_{\text{opp}}^{\,(n)},\vec{p}_{\text{self}}^{\,(n)}\odot(1-\vec{p}_{\text{opp}}^{\,(n)}),(1-\vec{p}_{\text{self}}^{\,(n)})\odot\vec{p}_{\text{opp}}^{\,(n)},(1-\vec{p}_{\text{self}}^{\,(n)})\odot(1-\vec{p}_{\text{opp}}^{\,(n)})\right)$
(2)

where $\odot$ denotes point-wise multiplication. Let $\vec{u}=(R,S,T,P)$
denote the vector of possible payoffs, $\vec{s}_{0}$ a one-hot vector
indicating the current state, and $\gamma_{\text{future}}<1$ a discount
factor. The average expected discounted payoff per round for
$\vec{p}_{\text{self}}$ is then given by

$\displaystyle U(\vec{p}_{\text{self}})$
$\displaystyle=\vec{s}_{0}^{T}\left(\sum_{k=1}^{\infty}\gamma_{\text{future}}^{k}T^{k}\right)\vec{u}/\left(\sum_{k=1}^{\infty}\gamma_{\text{future}}^{k}\right)$
$\displaystyle=(1-\gamma_{\text{future}})\vec{s}_{0}^{T}T\left(1-\gamma_{\text{future}}T\right)^{-1}\vec{u}\
.$ (3)

If the IPD has a fixed ending probability $p_{\text{end}}$ per step, we can
choose $\gamma_{\text{future}}=1-p_{\text{end}}$. However, in the rest of this
work, we will consider IPDs of fixed length, while assuming that the length
was not known to the participants beforehand. We will use
$\gamma_{\text{future}}=0.99$ in the rest of this work.

In order to find $\vec{p}_{\text{self}}$ which optimizes
$U(\vec{p}_{\text{self}})$, we use the Adam optimizer (Kingma and Ba,, 2014)
with a learning rate of 0.1 for 50 steps, using
$(\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2})$ as the starting point of
the optimization.

A conceptually similar strategy, IP0, was introduced by Lee et al., (2015).
IP0 optimizes payoffs in the stationary state against the memory-1 model of
the opponent, which is equivalent to Eq. (2.2) in the limit
$\gamma_{\text{future}}\rightarrow 1^{-}$.

###  2.3 CooperateISO: start cooperatively, then adapt

Finally, _CooperateISO_ combines the two sub-strategies introduced so far. It
starts by playing Longterm TFT, and then switches to playing ISO if it has
collected sufficient data about the opponent’s behavior (which gets expressed
in the memory-1 model) and responding optimally to the opponent model promises
higher payoffs than Longterm TFT has empirically achieved.

Note the trade-off involved in deciding when to start playing ISO. On the one
hand, playing Longterm TFT too long can be wasteful if a clearly better
response to the opponent exists. This includes opponents who do not retaliate
against defections, but also highly defective strategies, against which
Longterm TFT sub-optimally rewards each noise-induced cooperation. On the
other hand, switching to ISO is risky because it can destroy a cooperative
relationship without achieving anything better if our model of the opponent is
inaccurate. This can happen because we have not seen the opponent react in all
4 possible states sufficiently many times, because the opponent’s empirically
observed behavior is actually untypical for them because of noise events, or
because the opponent is simply not well described by a memory-1 model.

Note that in the presence of significant noise, maintaining cooperation with
unforgiving strategies like Grim Trigger (Banks and Sundaram,, 1990) is
hopeless irrespective of our behavior. This suggests that the threshold for
switching to ISO should be lower in the presence of high noise. However, for
simplicity we will use the same simple criterion for switching to ISO
irrespective of the noise level.

Formally, let $N_{c}$ be the number of rounds played so far (while playing
Longterm TFT), $U_{c}$ the average payoff per round achieved so far,
$\sigma_{c}$ the corresponding standard deviation, and $U_{a}$ the expected
average discounted payoff per round (from Eq. (2.2)) when playing ISO. We
switch from playing Longterm TFT to playing ISO if all of the following
conditions are met: $N_{c}\geq 10$ (require a minimum of data about the
opponent), $U_{a}-U_{c}>2\sigma_{c}/\sqrt{N_{c}}$ (the expected gain from
adapting is significant relative to the noise in the historical payoffs), and
$U_{a}-U_{c}>0.05(R-P)$ (the expected gain from adapting is non-negligible
compared to possible payoff differences).

###  2.4 Evaluation

The Axelrod library (Knight et al.,, 2016) contains (as of version 4.12.0) 239
strategies for the IPD from previous literature, including all strategies
mentioned so far and a zoo of different model types trained through
reinforcement learning (Harper et al.,, 2017). We use these strategies in
order to evaluate CooperateISO and its two sub-strategies, as well as compare
their performance with strategies which won previous tournaments or were
otherwise prominently discussed in previous literature. See Fig. 2 for this
evaluation.

Figure 2: Evaluation of the strategies introduced in Sec. 2 (CooperateISO and
its sub-strategies Longterm TFT and ISO) against the 239 strategies in the
Axelrod library, evaluated at noise levels of $0\%$, $1\%$, $5\%$, and $10\%$.
We use the conventional payoff values $T=5$, $R=3$, $P=1$, and $S=0$. Each IPD
lasts 400 steps, and each evaluated strategy plays 5 IPDs against each
opponent. Shaded regions show one standard error calculated as
$\sigma/\sqrt{5\cdot 400}$, where $\sigma$ is the sample standard deviation of
the average payoff per step.

Among the strategies _not_ introduced in the present work, we find
EvolvedLookerUp2_2_2 to be the strongest strategy in the noise-free case, and
DBS to be the strongest strategy in the presence of noise (with
$p_{\text{noise}}$ ranging from $1\%$ to $10\%$). This matches the findings of
Harper et al., (2017) (who used an older and smaller version of the Axelrod
library).

Longterm TFT by itself is not competitive with the best-performing strategies,
but outperforms the comparable strategies Generous TFT and Contrite TFT, in
particular in the presence of strong noise.

CooperateISO outperforms all previous strategies at all noise levels. ISO by
itself is clearly inferior to CooperateISO, in particular for zero or low
noise rates, where CooperateISO is much better able to maintain cooperation.
However, ISO outperforms its most similar strategy, DBS, despite doing simpler
opponent modelling. ISO’s peformance also approaches that of CooperateISO for
high noise levels, where even Longterm TFT has difficulties maintaining
cooperation, and so the difference between ISO and CooperateISO becomes
smaller.

##  3 Desiderata for IPD strategies: beyond tournaments

Tournaments provide valuable lessons about a strategy’s performance against a
wide range of opponents and served as the primary benchmark for evaluating
strategies in the IPD since Axelrod’s work in the 1980s (Axelrod,, 1984).
However, to the extend that we hope to draw wider-reaching lessons from
studying strategies in the IPD, tournament performance is an imperfect measure
of a strategy’s benefits. Firstly, a strategy’s performance in a tournament
always depends on the pool of strategies it competes against. Secondly, it
ignores aspects of the strategy that do not affect tournament performance but
will become relevant in other contexts.444 As a third reason, we note that the
ultimate goal in a tournament – achieving a good rank – is not a perfect proxy
for the goal in the IPD – achieving high payoff. For instance, the former
incentives submission of groups of strategies, in which all but one strategies
seek to push the lead strategy to a high rank (by always cooperating with it
and always defecting against all other participants). Such group strategies
indeed featured prominently in the 2004/2005 IPD tournaments (Kendall et al.,,
2007).  In particular, performance in a single tournament does not take into
account the effects of other parties analyzing our strategy and adjusting
their strategy to ours.

These considerations lead us to the following three (informal) desiderata,
which we motivate in more detail below. We stress that we are considering
scenarios in which we exclusively care about maximizing our strategy’s payoff.
In particular, we don’t care whether we achieve higher or lower payoffs than
the current opponent.555 Achieving higher payoffs than the opponent can be
crucial in population games, e.g. when seeking to resist an invading strategy
(Lee et al.,, 2015).

  1. 1.

The strategy should be _self-cooperating_ , i.e., achieve high rates of
cooperation when playing against a clone, including in the presence of
significant noise. In particular, when playing against a clone our strategy
should achieve an expected average payoff per step of $R-O(p_{\text{noise}})$
in the steady state. So a single noise event can only lead to $O(1)$
defections in self-play, and must not lead to a longer-lasting breakdown of
mutual cooperation.

  2. 2.

The strategy should be _cooperation-inducing_. That is, optimal (payoff-
maximizing) play against our strategy should lead to an expected average
payoff per step of $R-O(p_{\text{noise}})$ in the steady state for our
strategy.

  3. 3.

The strategy should be _adaptive_ (w.r.t. some set $\Omega$ of opponent
strategies). That is, our strategy should (be able to adapt to and thus)
achieve close-to-optimal expected payoffs against all strategies from some set
$\Omega$. We ware interested in sets $\Omega$ which consist of all
“sufficiently simple” strategies.

Let us briefly elaborate on why we chose these particular desiderata.

  * •

Not being self-cooperating limits our strategy’s payoff in any context in
which it has a high chance of facing a clone. This can happen, for example, if
a principal deploys multiple instances of the same strategy, or if other
players imitate our strategy.

  * •

If a strategy is not cooperation-inducing, it will do poorly as soon as an
opponent is able to create a decent model of it and react to that. TFT is the
simplest possible cooperation-inducing strategy, and the historically strong
performance of TFT and its generous variants lends pragmatic support to this
desideratum.

  * •

Note that there are strategies which satisfy a strictly stronger criterion
than being cooperation-inducing: there are “extortionate” zero-determinant
(ZD) strategies (Press and Dyson,, 2012) which have the property that optimal
(payoff-maximizing) play against them gives a payoff per step to them which is
higher than $R$. However, by construction in the IPD achieving a payoff higher
than $R$ per step implies that the opponent’s reward per step is lower than
$R$, so such extortionate strategies cannot be self-cooperating.

  * •

If a strategy is not adaptive w.r.t. some simple opponent, it is leaving “easy
money” on the table. For example, a weakness of TFT is that it always
cooperates against Cooperator, while a weakness of Pavlov is that it
cooperates on every second step against Defector. At the every least, we want
our strategy to be adaptive w.r.t. all memory-0 strategies. That is, our
strategy should always defect against all strategies which cooperate with
constant probability irrespective of past actions and thus provide no
incentive to cooperate. DBS is adaptive w.r.t. all memory-1 strategies666 It
_tries_ to learn a memory-1 model of the opponent and react optimally to that.
We leave aside the more difficult question of how reliably it achieves that in
practice and in the presence of noise. , and its strong performance in
tournaments lends pragmatic support to this desideratum.

  * •

Note that we do not add _evolutionary stability_ as a desideratum – there are
no evolutionary stable strategies in the IPD (Selten and Hammerstein,, 1984;
Boyd and Lorberbaum,, 1987; Farrell and Ware,, 1989; Lorberbaum,, 1994).
Instead, if the set of strategies is not restricted, evolutionary dynamics
will move populations between a variety of Nash equilibria with different
levels of cooperation (García and van Veelen,, 2018).

Similar desiderata have been proposed by Neill, (2001). They search for
strategies which are “self-cooperating” (able to achieve mutual cooperation
with their clone), “C-exploiting” (able to exploit unconditional cooperators),
and “D-unexploitable” (able to resist exploitation by defectors). Using our
terminology, the latter two correspond to being adaptive w.r.t.
$\Omega=\\{\text{Cooperator},\text{Defector}\\}$, the weakest non-trivial form
of adaptiveness. We also note that being “D-unexploitable” is a corollary of
being cooperation-inducing.

Similar desiderata have also been discussed in the context of population
dynamics. IP0 (Lee et al.,, 2015) is adaptive w.r.t. all memory-1 strategies
(as discussed in Sec. 2.2), and self-cooperating by virtue of a noise-tolerant
handshake mechanism. These properties make IP0 both uninvadable and a strong
invader when evaluated against memory-1 strategies. Knight et al., (2018)
found that strategies which are best at invading populations of other
strategies tend to be strategies trained to do well in tournaments; while the
best resistors of invasion invoke handshake mechanisms to cooperate with each
other, but not with invaders.

Table 1 assesses fulfillment of the proposed desiderata by prominent
strategies from previous literature as well as the strategies introduced in
the present work. Fig. 3 empirically evaluates self-cooperativeness at
$p_{\text{noise}}=5\%$. CooperateISORevert1 and CooperateISORevert2 are
introduced in the following section.

| Self-cooperating | Cooperation-inducing | Adaptive
---|---|---|---
Pavlov | ✓ | ✗777In the limit of vanishing noise, always cooperating gets a payoff per step of $R$ against Pavlov while always defecting gets $(T+P)/2$. With the conventional IPD payoffs, both of these are equal. With finite noise and the conventional IPD payoffs, always defecting becomes preferable. | ✗
TFT, Omega TFT | ✗ | ✓ | ✗
Generous/Contrite TFT | ✓ | ✓ | ✗
DBS | ✗ | ✗888The way DBS updates its opponent model makes it non-cooperation-inducing. DBS creates a deterministic rule (of the form “the opponent always cooperates/defects in this state”) once the opponent plays the same action 3 times consecutively in the given state. It only abandons such a rule after 3 consecutive violations; 2 consecutive violations are forgotten if the next action follows the rule. This makes it easy to exploit DBS: play in accordance with TFT if that is necessary to ensure that DBS’ model of the opponent is TFT (and hence DBS cooperates, since that is the best response against TFT); otherwise defect (and rely on DBS ignoring up to 2 consecutive violations of its opponent model). Empirically, we find that such a simple exploitative strategy can achieve rewards per step much higher than $R$, leading to rewards much lower than $R$ for DBS. | ✓(1)
IP0 | ✓ | ✗ | ✓(1)
Longterm TFT | ✓ | ✓ | ✗
ISO, CooperateISO | ✗ | ✗ | ✓(1)
CooperateISORevert1 | ✓ | ✗ | ✓(1)
CooperateISORevert2 | ✓ | ✓ | ✓(2)
Table 1: Fulfillment of the three proposed desiderata by prominent or newly
introduced strategies. The second half of the table describes strategies
introduced in the present work. For adaptiveness, ✓(1) denotes being adaptive
w.r.t. all memory-1 strategies, and ✓(2) denotes being adaptive w.r.t. all
_non-extortionate_ memory-1 strategies.

Figure 3: Average payoff per step achieved by various strategies in self-play
(playing against a clone) in the IPD, using $p_{\text{noise}}=5\%$. The top
figure shows strategies from previous literature, the bottom figure shows
strategies introduced in the present work. Each payoff is averaged over 100
self-play games. The numbers in parentheses after each strategy name show the
average payoff per step achieved by that strategy (averaged over 1000 steps
per game and 100 games). Dashed lines show three simple benchmarks: always
cooperate, cooperate/defect randomly with $50\%$ probability, always defect,
while including the effects of noise.

If we accept these desiderata, a natural next question is whether it is
possible for a strategy to fulfill all of them; or, more precisely, what the
“maximal” set $\Omega$ is such that it is possible for a strategy to fulfill
all of them. We note that it is not possible for a strategy to be both
cooperation-inducing and adaptive w.r.t. all memory-1 strategies. This follows
from the existence of extortionate ZD strategies (Press and Dyson,, 2012),
which are memory-1 strategies and have the property that optimally responding
to them (in the sense of maximizing payoffs) yields an average payoff higher
than $R$ per step for them. So if a strategy is adaptive w.r.t. all memory-1
strategies, optimal play against that strategy can achieve an average payoff
higher than $R$ per step, implying less than $R$ per step for the strategy.
Let us define a strategy as _extortionate_ if the payoff-maximizing response
to the strategy gives average payoffs higher than $R$ to the strategy (which
implies less than $R$ for the opponent). The best we can hope for is thus a
strategy which is self-cooperating, cooperation-inducing, and adaptive to all
non-extortionate memory-1 strategies. The following section introduces such a
strategy (to the best of our knowledge, the first such strategy).

##  4 Cooperate, adapt, and revert: fulfilling all desiderata

In order to create a strategy which is self-cooperating, cooperation-inducing,
and adaptive to all non-extortionate memory-1 strategies, we start with
CooperateISO and add two rules to it for reverting from the state in which it
plays ISO to the state in which it plays Longterm TFT. The basic idea is
illustrated in Fig. 1.

###  4.1 Becoming self-cooperating: revert after poor performance

CooperateISO is adaptive, but neither self-cooperating nor cooperation-
inducing. When playing against itself, it will initially be highly forgiving
as long as it is in the Longterm TFT state. Optimal play against a forgiving
opponent is to defect. Both clones will thus switch to playing ISO and always
defecting. They will then update their model of the opponent to being highly
defective, against which the optimal response is to defect further. (DBS goes
through a similar sequence of updates when playing against a clone.)

In order to break this cycle, we add a rule to CooperateISO to revert from the
ISO state to the Longterm TFT state if the former empirically performs worse
than the latter. This ensures self-cooperativeness, and might also seem like a
prudent choice when playing against other opponents. We call the resulting
strategy _CooperateISORevert1_. In order to decide whether ISO’s payoffs are
lower than the ones of Longterm TFT, we use a standard significance test.
Formally, let $N_{c}$ the number of rounds for which we played Longterm TFT,
$U_{c}$ the average payoff per step, and $\sigma_{c}$ the corresponding
standard deviation; and let $N_{a}$, $U_{a}$, and $\sigma_{a}$ be the
analogous values while playing ISO. We revert from playing ISO to playing
Longterm TFT if $N_{a}\geq 10$ and

$\displaystyle\frac{U_{c}-U_{a}}{\sqrt{\frac{\sigma_{c}^{2}}{N_{c}}+\frac{\sigma_{a}^{2}}{N_{a}}}}>2\
.$ (4)

The effect of this additional rule can be seen in the bottom plot in Fig. 3.
When CooperateISO plays against a clone, payoffs degrade to far below the
level of mutual cooperation. For CooperateISORevert1 (and CooperateISORevert2,
which will be introduced in the following subsection), payoffs degrade
similarly for the first few dozen steps. They subsequently recover when the
strategy is likely to have reverted to playing Longterm TFT.

###  4.2 Becoming cooperation-inducing: revert if extorted

CooperateISORevert1 is still not cooperation-inducing; an extortionate ZD
strategy can still get rewards per step higher than $R$, if ISO responds
optimally to it.999 Extortionate strategies face difficulties in the presence
of noise, however; see Hao et al., (2015) for a discussion.  In order to
create a strategy which is cooperation-inducing, we add a second rule for
reverting from ISO to Longterm TFT: revert to Longterm TFT if ISO is being
extorted, which we define as the opponent achieving payoffs per step which are
higher than $R$. We call the resulting strategy _CooperateISORevert2_.
Formally, let let $N_{o}$, $U_{o}$, and $\sigma_{o}$ as above describe the
payoffs which the opponent achieves while playing against ISO. We revert to
playing Longterm TFT if $N_{o}\geq 10$ and $U_{o}-2\sigma_{o}/\sqrt{N_{o}}>R$.

By construction, it is not possible to achieve payoffs larger than $R$ in the
steady state against CooperateISORevert2 while CooperateISORevert2 is in the
ISO state. TFT is the least cooperative strategy which ensures Longterm TFT’s
cooperation in the steady state, and so is the best response to Longterm TFT.
The average payoff per step which TFT achieves against Longterm TFT is
$R+(S+2T-3R)p_{\text{noise}}+O(p_{\text{noise}}^{2})$, since the two leading-
order effects of noise are Longterm TFT accidentally defecting and TFT
retaliating, and TFT accidentally defecting with Longterm TFT _not_
retaliating. With the conventional payoffs, $S+2T-3R>0$. TFT thus also
achieves payoffs of $R+(S+2T-3R)p_{\text{noise}}+O(p_{\text{noise}}^{2})\geq
R$ against CooperateISORevert2, making it the best response to
CooperateISORevert2. This makes CooperateISORevert2 cooperation-inducing.

Figure 4: Evaluation of all strategies introduced in the present work against
the 239 strategies in the Axelrod library, evaluated at noise levels of $0\%$,
$1\%$, $5\%$, and $10\%$. “Best previous” shows the performance of
EvolvedLookerUp2_2_2 for $p_{\text{noise}}=0$ and DBS for
$p_{\text{noise}}>0$. Each IPD lasts 400 steps, and each evaluated strategy
plays 5 IPDs against each opponent. Shaded regions show one standard error
calculated as $\sigma/\sqrt{5\cdot 400}$, where $\sigma$ is the sample
standard deviation of the average payoff per step.

Fig. 4 evaluates all strategies developed in the present work against the
strategies in the Axelrod library. We also compare all of them to the
previously existing strategies which, to the best of our knowledge, show the
strongest performance: EvolvedLookerUp2_2_2 for $p_{\text{noise}}=0$ and DBS
for $p_{\text{noise}}>0$. Ensuring that the strategy is both self-cooperating
and cooperation-inducing, i.e., playing CooperateISORevert2 instead of
CooperateISO, has a small cost in performance against this pool of opponents.
However, CooperateISORevert2 still outperforms the best-performing previously
existing strategies at all noise levels.

##  5 Discussion

While even some bacteria show TFT-like behavior (Smith et al.,, 2020), the
strategies introduced in the present work are arguably too complex to evolve
in biological systems. The lessons from this work are thus more relevant in
the context of humans, human organizations, and human-designed systems which
interact in iterated social dilemmas. In such contexts, it is plausible that
opponents imitate or analyze our strategy and adapt theirs to it. This makes
it desirable that our strategy be both self-cooperating and cooperation-
inducing.

When facing a diverse set of opponents, CooperateISO shows the value of
enabling robust mutual cooperation, while also being ready to adapt oneself to
the opponent if it does not incentivize cooperative behavior. We have shown
that such a strategy can be updated to one which is both self-cooperating and
cooperation-inducing, without incurring a large loss in performance even in
the more narrow context of a tournament against a fixed set of opponents.

## Acknowledgements

The author thanks Marc Harper for valuable comments on this manuscript.

## References

  * Au and Nau, (2007) Au, T.-C. and Nau, D. (2007).  Is it accidental or intentional? a symbolic approach to the noisy iterated prisoner’s dilemma.  The Iterated Prisoners Dilemma, 20:231–262.
  * Axelrod, (1984) Axelrod, R. (1984).  The Evolution of Cooperation.  Basic Books.
  * Axelrod and Hamilton, (1981) Axelrod, R. and Hamilton, W. D. (1981).  The evolution of cooperation.  Science, 211(4489):1390–1396.
  * Banks and Sundaram, (1990) Banks, J. S. and Sundaram, R. K. (1990).  Repeated games, finite automata, and complexity.  Games and Economic Behavior, 2(2):97–117.
  * Boyd and Lorberbaum, (1987) Boyd, R. and Lorberbaum, J. (1987).  No pure strategy is evolutionarily stable in the repeated prisoner’s dilemma game.  Nature, 327:58–59.
  * Farrell and Ware, (1989) Farrell, J. and Ware, R. (1989).  Evolutionary stability in the repeated prisoner’s dilemma.  Theoretical Population Biology, 36(2):161–166.
  * García and van Veelen, (2018) García, J. and van Veelen, M. (2018).  No strategy can win in the repeated prisoner’s dilemma: Linking game theory and computer simulations.  Frontiers in Robotics and AI, 5:102.
  * Hao et al., (2015) Hao, D., Rong, Z., and Zhou, T. (2015).  Extortion under uncertainty: Zero-determinant strategies in noisy games.  Physical Review E, 91(5):052803.
  * Harper et al., (2017) Harper, M., Knight, V., Jones, M., Koutsovoulos, G., Glynatsi, N. E., and Campbell, O. (2017).  Reinforcement learning produces dominant strategies for the iterated prisoner’s dilemma.  PloS one, 12(12):e0188046.
  * Kendall et al., (2007) Kendall, G., Yao, X., and Chong, S. Y. (2007).  The iterated prisoners’ dilemma: 20 years on, volume 4.  World Scientific.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Knight et al., (2016) Knight, V., Campbell, O., Harper, M., Langner, K., Campbell, J., Campbell, T., Carney, A., Chorley, M., Davidson-Pilon, C., Glass, K., et al. (2016).  An open reproducible framework for the study of the iterated prisoner’s dilemma.  arXiv preprint arXiv:1604.00896.
  * Knight et al., (2018) Knight, V., Harper, M., Glynatsi, N. E., and Campbell, O. (2018).  Evolution reinforces cooperation with the emergence of self-recognition mechanisms: An empirical study of strategies in the moran process for the iterated prisoner’s dilemma.  PloS one, 13(10):e0204981.
  * Lee et al., (2015) Lee, C., Harper, M., and Fryer, D. (2015).  The art of war: Beyond memory-one strategies in population games.  PloS one, 10(3):e0120625.
  * Lorberbaum, (1994) Lorberbaum, J. (1994).  No strategy is evolutionarily stable in the repeated prisoner’s dilemma.  J Theor Biol., 168(2):117–30.
  * Molander, (1985) Molander, P. (1985).  The optimal level of generosity in a selfish, uncertain environment.  Journal of Conflict Resolution, 29(4):611–618.
  * Neill, (2001) Neill, D. B. (2001).  Optimality under noise: higher memory strategies for the alternating prisoner’s dilemma.  Journal of Theoretical Biology, 211(2):159–180.
  * Nowak and Sigmund, (1993) Nowak, M. and Sigmund, K. (1993).  A strategy of win-stay, lose-shift that outperforms tit-for-tat in the prisoner’s dilemma game.  Nature, 364(6432):56–58.
  * Nowak, (2006) Nowak, M. A. (2006).  Five rules for the evolution of cooperation.  science, 314(5805):1560–1563.
  * Press and Dyson, (2012) Press, W. H. and Dyson, F. J. (2012).  Iterated prisoner’s dilemma contains strategies that dominate any evolutionary opponent.  Proceedings of the National Academy of Sciences, 109(26):10409–10413.
  * Selten and Hammerstein, (1984) Selten, R. and Hammerstein, P. (1984).  Gaps in harley’s argument on evolutionarily stable learning rules and in the logic of “tit for tat”.  Behav. Brain Sci., 7:115–116.
  * Slany et al., (2007) Slany, W., Kienreich, W., et al. (2007).  On some winning strategies for the iterated prisoner’s dilemma, or, mr. nice guy and the cosa nostra.  The iterated prisoners’ dilemma, 20:171.
  * Smith et al., (2020) Smith, W. P., Brodmann, M., Unterweger, D., Davit, Y., Comstock, L. E., Basler, M., and Foster, K. R. (2020).  The evolution of tit-for-tat in bacteria via the type vi secretion system.  Nature communications, 11(1):1–11.
  * Wu and Axelrod, (1995) Wu, J. and Axelrod, R. (1995).  How to cope with noise in the iterated prisoner’s dilemma.  Journal of Conflict resolution, 39(1):183–189.


# CNN-tuned spatial filters for P- and S-wave decomposition and applications
in elastic imaging

Wenlong Wang3, Jianwei Ma4, George A. McMechan6
3Center of Geophysics, Department of Mathematics and Artificial Intelligence
Laboratory, Harbin Institute of Technology, Harbin, China 150001
(e-mail:wenlong.wang@hit.edu.cn).  4School of Earth and Space Sciences, Peking
University, Beijing, China, 100871.  6Center for Lithospheric Studies, The
University of Texas at Dallas, 800 W Campbell Road (ROC21), Richardson, TX,
USA 75080.

###### Abstract

P- and S-wave decomposition is essential for imaging multi-component seismic
data in elastic media. A data-driven workflow is proposed to obtain a set of
spatial filters that are highly accurate and artifact-free in decomposing the
P- and S-waves in 2D isotropic elastic wavefields. The filters are formulated
initially by inverse Fourier transforms of the wavenumber-domain operators,
and then are tuned in a convolutional neural network to improve accuracy using
synthetic snapshots. Spatial filters are flexible for decomposing P-and
S-waves at any time step without performing Fourier transforms over the entire
wavefield snapshots, and thus are suitable for target-oriented imaging.
Snapshots from synthetic data show that the network-tuned spatial filters can
decompose P- and S-waves with improved decomposition accuracy compared with
other space-domain PS decomposition methods. Elastic reverse-time migration
using P- and S-waves decomposed from the proposed algorithm shows reduced
artifacts in the presence of a high velocity contrast.

###### Index Terms:

PS decomposition, Convolutional neural network, Elastic imaging

##  I Introduction

P- and S-waves coexist in seismic waves propagating through the solid earth.
The causes, polarization directions, and propagation velocities, of P- and
S-waves are different, and thus they carry different information from the
subsurface at a variety of scales. For example, in seismic imaging for oil and
gas exploration, P- and S-wave reflections can be imaged independently. In the
presence of gas clouds and fractures, the P-waves are strongly attenuated,
while S-waves are less affected, and S-wave images provide valuable structural
information from reservoirs [1, 2]. In computational seismology, most imaging
methods for earthquake and microseismic events involve picking the P- and
S-wave arrival times. Artman _et al._ [3] and Yang and Zhu [4] propose time-
reversal imaging of earthquake ruptures and hydraulic fractures by separating
and propagating the recorded P- and S-waves back into the earth, and applying
the elastic imaging conditions. All of these applications require separating
the P- and S-waves accurately and efficiently.

One approach to separate and utilize both P- and S-wave data is the ray-based
method. Examples include Kirchhoff migration [5, 6, 7], in which the PS
separation is implicit by using P- and S-velocity models separately for ray-
tracing to compute the image times. Ray-based methods are efficient, but their
high frequency assumption makes ray theory unable to fully describe wave
phenomena, especially for complicated models [8]. Another category is wave-
equation based solutions [9, 10, 11], which reconstruct the full elastic
receiver wavefields from boundary conditions [12], and separates the P- and
S-waves before, or as a part of, application of the image condition [13].

The divergence and curl operators, which are based on Helmholtz decomposition
[14], are widely applied in PS separations. However, the phases and amplitudes
of the separated P- and S-waves are distorted by the divergence and curl
operators, causing problems for true amplitude imaging and subsequent
interpretations [15]. To solve this problem, PS wavefield decomposition is
proposed, which preserves the vector components of the input elastic wavefield
[16, 17, 18].

There are two groups of algorithms for PS decomposition depending on the
domain of calculation. In the wavenumber domain, elastic wavefields are
projected to the P- and S-polarization directions [17]; Zhu [19] decomposes
the wavefields by solving a Poisson’s equation. Both of these methods involve
Fourier transforms. However, it is expensive, or even prohibitive, to perform
Fourier transforms on large (3D) datasets. Xiao and Leaney [20] introduce an
auxiliary stress wavefield into the elastodynamic wave equations to obtain
decoupled wavefields, but the decomposition is embedded in the elastic
extrapolations (decoupled propagation), and artifacts are likely to be
generated at large velocity contrasts in the model [18]. Yan and Sava [13]
propose spatial filters for separating P- and S-waves in transversely
isotropic media with a vertical axis (VTI). The spatial filters don’t require
Fourier transforms, but the accuracy is limited by the truncation of the
spatial filters.

PS separation and decomposition can also be achieved with deep learning
algorithms. Wei et al. [21] reconstruct P-waves in VSP data using conditional
generative adversarial network. Wei et al. [22] extend the algorithm to
separate P- and S-waves in multi-component VSP data using multi-task learning.
However, these methods need a lot of training samples from a variety of
velocity models to compensate the lack of physical information. The spatial
filters are convolutional operators, which are widely used in the computer
vision community [i.e. convolutional neural networks (CNNs)] [23]. Thus, CNNs
are suited for separating P- and S-wave snapshots. Kaur _et al._ [24] apply a
generative adversarial network (GAN) to match the predicted P-waves with the
decomposed ones from low-rank approximation in anisotropic media; the network
needs to be re-trained to apply to different models. We propose to formulate
the PS decomposition as a single-layer CNN with the coupled wavefield
snapshots as inputs and the decomposed P-wave snapshots as outputs. The
spatial filters are tuned to optimize PS decomposition accuracy by training in
a training set. Different from the existing methods, our network only needs to
be trained once on a complex model and can be applied to any models, because
the physics in PS decomposition are embedded in the network structure. The
trained network can decompose elastic wavefields in different models without
re-training.

This paper is organized as follows, we first illustrate the wavenumber domain
operators for P- and S-wave decomposition, and their relation to spatial
filters. Then we describe a CNN structure to tune the spatial filters.
Wavefield snapshot decomposition and elastic reverse time migration are
performed to document the improvements in accuracy by comparing the
decomposition results with those from two other space-domain PS decomposition
methods. To limit the scope, we focus on 2D isotropic media, although PS
decompositions in anisotropic and anelastic media are also possible [25, 26,
27].

##  II PS decomposition

Different from the commonly used divergence and curl operators [14] which
distort the amplitudes and phases during wavefield separation, PS
decomposition preserves all the components of the elastic vector wavefield. In
the wavenumber domain, the explicit equations to obtain 2D decomposed P-waves
are [17]

$\tilde{U}^{P}_{x}=K^{2}_{x}\tilde{U}_{x}+K_{x}K_{z}\tilde{U}_{z},$ (1)

and

$\tilde{U}^{P}_{z}=K^{2}_{z}\tilde{U}_{z}+K_{z}K_{x}\tilde{U}_{x},$ (2)

where $K_{x}$ and $K_{z}$ are the components of the normalized wavenumber
vector $\bm{K}$, which represents the P-wave polarization direction in
isotropic media. $U_{x}$ and $U_{z}$ are the space-domain elastic wavefield
components with coupled P- and S-waves, and a tilde over the wavefield
indicates the wavenumber-domain representation. The S-waves can be decomposed
in a similar way, or by subtracting the P-waves from the coupled wavefields,
component-by-component-by

$\tilde{U}^{S}_{x}=\tilde{U}_{x}-\tilde{U}^{P}_{x},$ (3)

and

$\tilde{U}^{S}_{z}=\tilde{U}_{z}-\tilde{U}^{P}_{z},$ (4)

in either the space or the wavenumber domain. Thus only the P-wave
decomposition is discussed and demonstrated in the following parts of this
paper.

##  III Spatial filters for P- and S-wave decomposition

In the 2D equations 1 and 2, three wavenumber operators $K^{2}_{x}$,
$K^{2}_{z}$ and $K_{x}K_{z}$ are involved; the $K^{2}_{z}$ is the transpose of
$K^{2}_{x}$ when the spatial intervals along $x$ and $z$ directions are the
same. $K^{2}_{x}$ and $K_{x}K_{z}$ are plotted in Figures 1a-1b, respectively.
The wavenumber domain operators can be transformed into the space-domain by
inverse Fourier transforms to formulate the spatial filters [13]. The
equivalent expressions of equations 1 and 2 in the space-domain are

$U^{P}_{x}=L_{x}[U_{x}]+L_{xz}[U_{z}],$ (5)

and

$U^{P}_{z}=L_{z}[U_{z}]+L_{xz}[U_{x}],$ (6)

where $L_{x}$, $L_{z}$ and $L_{xz}$ are filters that are designed to decompose
the elastic wavefields. The square brackets represent spatial filtering of the
wavefields.

Figure 1: The wavenumber domain operators (a) $K^{2}_{x}$ and (b) $K_{x}K_{z}$
for isotropic P- and S-wave decompositions.

Different from the anisotropic media, the polarizations of P- and S-waves in
isotropic wavefield have much simpler relations with the wavenumbers, they do
not depend on the model parameters ($V_{P}$, $V_{S}$ or $\rho$). Thus the
spatial filters are stationary and universal for different models. The spatial
filters are trained only once using synthetic labelled data pairs, and the
trained filters can be applied to decompose P- and S-waves in any isotropic
wavefields. Full representations of the wavenumber domain operators are
infinitely large 2D patches in the space-domain which make them infeasible to
be applied. Yan and Sava [13] propose to truncate the filters to balance the
computation efficiency and the separation accuracy. A set of truncated $L_{x}$
and $L_{xz}$ filters of size 15 $\times$ 15 are shown in Figures 2a-2b
($L_{z}$ is the transpose of $L_{x}$, and thus is not plotted).

Figure 2: The (a) $L_{x}$ and (b) $L_{xz}$ spatial filters, which are
transformed from the wavenumber domain operators $K^{2}_{x}$ and $K_{x}K_{z}$
in Figure 1. $L_{z}$ is the transpose of $L_{x}$.

##  IV Tuning the spatial filters with CNN

Truncations of the filters reduce the accuracy and completeness of the
decomposition. To improve the accuracy without increasing the size of the
filters, we propose to tune the spatial filters $L_{x}$ and $L_{xz}$ by
training in a CNN. Following equations 5 and 6, a CNN is constructed to obtain
the P-waves $\bm{U}^{P}$ from the coupled wavefield snapshot $\bm{U}$ of the
form

$\bm{U}^{P}=\mathrm{CNN}(\bm{U}).$ (7)

The two components ($U_{x}$ and $U_{z}$) of the elastic wavefield snapshot
$\bm{U}$ are input to the network as two channels, and the output is expected
to be the decomposed 2-component P-wave snapshot $\bm{U}^{P}$.

The CNN for PS decomposition is shown in Figure 3. Two spatial filters
($L_{x}$ and $L_{xz}$) are to be trained. Most neural networks in deep
learning are black-boxes and cannot be explained. However, the CNN which we
use contains only one CNN layer, which can be simply interpreted as spatial
filters. The physical meaning of our CNN is a spatial representation of a
Fourier domain operator, which is similar to the work of Yan and Sava [13].
Thus our network does not contain pooling or non-linear activation functions
because of the well understood physics that is involved. Another difference is
that, instead of initializing the spatial filters with random numbers, they
are obtained by transforming the wavenumber domain operators in Figure 2a and
2b, and are set as initial values before training.

The network is trained by minimizing the loss function

$L=\sum_{\begin{subarray}{c}i=1,N\end{subarray}}|\bm{U}^{P}_{i}-\hat{\bm{U}}^{P}_{i}|^{2},$
(8)

where $\bm{U}^{P}$ and $\hat{\bm{U}}^{P}$ are the predicted and ground truth
P-waves, respectively. $|\cdot|$ calculates the modulus arguments of vector
difference. $N$ is the number of grid points in a wavefield snapshot. The
filters are tuned in a CNN training process. We expect the training to reduce
the errors that are associated with the filter truncations.

We generate a training set with 800 samples, which contains coupled and
decomposed P-wave snapshot pairs. The coupled wavefield snapshots are
generated by simulating wave propagation in a portion of the Sigsbee model
[28] (Figure 4) using an eighth-order in space, second-order in time, stress-
particle-velocity, staggered-grid, finite-difference solution [29]. An
explosive source with a 10 Hz Ricker wavelet is placed at (x, z) = (1.28,
0.02) km. The spatial intervals in both x and z directions are 10 m, and the
time interval is 1 ms. The 800 2-component snapshots used are captured at
random time steps from T = 0.0 s to T = 1.4 s. Although reflections from model
edges can also be separated into P- and S-waves, we use Fourier transforms to
generate the ground truth (decomposed P- and S-wave snapthots), so they may
have wrap-around errors at the model edges. Thus, convolutional perfectly
matched layer (CPML) absorbing boundary conditions [30] are used on all four
boundaries to reduce unwanted reflections.

The training doesn’t need a massive training set, because (1) approximate
spatial filters (Figures 2a and 2b) are set as initial filters which are close
to the optimized filters; (2) the CNN structure in Figure 3 imposes a physical
contraint on the training; (3) only 2 filters ($L_{x}$ and $L_{xz}$) need to
be tuned. After training, the CNN-tuned spatial filters can be applied to any
wavefield snapshot that has the same spatial intervals as the training set.

The ground truth P-waves are the decomposition results using the wavenumber
decomposition algorithm (equations 1 and 2), which is, by far, the most
accurate algorithm for P- and S-wave decomposition [17]. Samples of the
training set are plotted in Figure 5. The training set is sufficiently
diverse, as the Sigsbee model contains a wide range of structures and
velocities.

Figure 3:  The CNN architecture with respect to equations 5 and 6 for tuning
the spatial filters. The arrows indicate the data flow direction, and the
$L_{x}$, $L_{z}$ and $L_{xz}$ are convolutional operators that operate on the
input wavefields. The $L_{z}$ can be obtained by transposing $L_{x}$, and thus
they share the same filter values during training.  Figure 4:  The (a) Vp, (b)
Vs and (c) density of the Sigsbee model for generating wavefield snapshots for
training.  Figure 5:  Four representative snapshot pairs in the training set.
(a) The coupled X and Z particle velocity components and (b) the corresponding
(ground truth) decomposed P-wave x and z components obtained using equations
1-2.

The size of the spatial filter should be at least one wavelength to maximize
the resolution, but large spatial filters increase the computation cost. Three
sets of spatial filters are trained with dimension 9 $\times$ 9, 15 $\times$
15 and 21 $\times$ 21\. The networks are trained using an Adam optimizer [31]
in the Pytorch (www.pytorch.org) platform. Figure 6 shows that losses decrease
as the training proceeds, and that the decomposition accuracy of all these
sizes of filters is improved during training. As expected, spatial filters
with larger sizes can reach lower losses. Big improvements can be made when
the filters expand from 9 $\times$ 9 to 15 $\times$ 15, but improvements are
less significant for filter size 21 $\times$ 21 (Figure 6). In the following
tests, only the updated $L_{x}$ and $L_{xz}$ filters of the dimension 15
$\times$ 15 are used, which are plotted in Figures 7a and 7b, respectively.
Figures 7c and d show the differences between the filters before and after the
tuning.

Figure 6:  The loss evolutions of different filter sizes during traning.
Figure 7:  The (a) $L_{x}$ and (b) $L_{xz}$ spatial filters after tuning with
CNN. (c) and (d) show the difference between the $L_{x}$ and $L_{xz}$ filters
before and after tuning, respectively.

##  V Synthetic tests and comparisons

We perform tests on synthetic data using the CNN-tuned spatial filters (CNN-
SF), and the results are compared with the decomposed P-waves from untuned
spatial filters (SF), and those from decoupled propagation (DP) [16]. DP
involves solving an auxiliary equation along with the elastodynamic wave
equations, and thus has to be applied at each time step from the beginning of
the elastic wavefield extrapolation. Note that CNN-SF is tuned only from the
Sigsbee model in the previous section, and can be directly applied to the
following (different) models without re-training.

###  V-A Two-layer example

The first test is performed on an isotropic two-layer model with $V_{P}$ = 3
km/s, $V_{S}$ =2.1 km/s, and $\rho$ = 2.2 $\mathrm{g/cm^{3}}$ for the upper
layer, and $V_{P}$ = 4 km/s, $V_{S}$ =2.4 km/s, and $\rho$ = 2.4
$\mathrm{g/cm^{3}}$ for the lower layer. The model has 10 m grid spacing in
both the x- and z-directions. An explosive source with a 10 Hz Ricker wavelet
is placed at (x, z) = (1.28, 0.9) km. A snapshot of the PS coupled x and z
particle velocity components at t = 0.42 s is shown in Figure 8a, and the
ground truth P-wave using equations 1 and 2 are shown in Figure 8b as
benchmarks.

Figure 8:  The (a) coupled x and z wavefield components and (b) decomposed
P-wave x and z components (the ground truth) in the two-layer model at t =
0.42 s.

The decomposed P-waves using DP, SF and CNN-SF methods are shown in Figures
9a-9c, respectively. To measure the decomposition accuracy, we subtract the
decomposition results ($\bm{U}^{P}_{i}$) in Figures 9a-9c from the ground
truth ($\hat{\bm{U}}^{P}_{i}$) P-wave (Figure 8b), and the residuals are shown
in Figures 9d-9f. DP has high accuracy at homogeneous parts of the medium, but
it generates serious artifacts along the interface between layers, which are
caused by the added auxiliary equation [18], indicating that the DP method can
be applied only to smooth models. The spatial filter algorithms SF and CNN-SF
don’t have this problem. The residuals in Figure 9f are smaller than in Figure
9e, indicating that the CNN-SF achieves higher accuracy than the SF. We
quantify the accuracy by

$A=1-\frac{\sum_{\begin{subarray}{c}i=1,N\end{subarray}}|\bm{U}^{P}_{i}-\hat{\bm{U}}^{P}_{i}|^{2}}{\sum_{\begin{subarray}{c}i=1,N\end{subarray}}|\bm{U}^{P}_{i}|^{2}}.$
(9)

The calculated accuracies of DP, SF and CNN-SF in this example are 92.2%,
90.1%, and 98.6%, respectively.

Figure 9:  The PS decomposition results using (a) DP, (b) SF and (c) CNN-SF in
the two-layer model. The decomposition accuracy can be measured by subtracting
their respective decomposition results from the ground truth (Figure 8b); the
residuals are shown in (d), (e) and (f), respectively.

The computation times of DP, SF, and CNN-SF are 0.01, 0.25 and 0.25 seconds,
respectively, for each decomposition of the above model. DP has the best
efficiency, but it involves solving an auxiliary equation along with solving
the wave equations, thus decomposition with DP is not independent from
wavefield extrapolation. CNN-SF and SF algorithms, on the other hand, don’t
require extrapolations, and they have the same filter size, and thus share the
same computation time.

###  V-B Marmousi-2 model example

The second test is performed on a portion of the Marmousi-2 model [32] (Figure
10). The model has 10 m grid spacing in both x- and z-directions. An explosive
source with a 15 Hz Ricker wavelet is placed at (x, z) = (1.28, 0.04) km. A
snapshot of the x and z particle velocity components at t = 0.21 s is shown in
Figure 11.

Figure 10:  The (a) P-wave velocity, (b) S-wave velocity and (c) density of
the Marmousi-2 model.  Figure 11:  The (a) coupled wavefield snapshot X and Z
components and (b) decomposed P-wave snapshot X and Z components (ground
truth) in the Marmousi-2 model at t = 0.21 s.

We decompose the representative coupled wavefields in Figure 11 using the DP,
SF and CNN-SF methods, respectively; the decomposed P-waves are shown in
Figures 12a-12c, and their residuals are in Figures 12d-f, respectively.
Results of DP contain artifacts along the velocity interfaces. The
decomposition results from SF and CNN-SF methods don’t have such artifacts,
and the decomposition accuracy of CNN-SF is better than SF.

Figure 12:  The PS decomposition results using (a) DP, (b) SF and (c) CNN-SF
in the Marmousi-2 model. The decomposition residuals are shown in (d), (e) and
(f), respectively.

To have a comprehensive analysis of the decomposition accuracy, 575 snapshots
of coupled wavefields at different time steps are captured from the wavefield
simulation in the Marmousi-2 model, and we use DP, SF and CNN-SF to generate
the decomposed P-waves separately from all the snapshots. The accuracy of the
three algorithms with increasing time steps are plotted in Figure 13. As the
wave propagates through the model, more converted S-waves are generated, and
the accuracies of all three algorithms decrease. Because the Marmousi-2 model
has more velocity interfaces than the two-layer model, the accuracy of DP
decreases significantly as the wave propagates deeper into the model. CNN-SF
has the best accuracy among the three algorithms at all time steps.

Figure 13:  Decomposition accuracy of DP, SF and CNN-SF algorithms of
wavefield snapshots at different extrapolation time steps.

###  V-C Application to elastic reverse-time migration

In elastic reverse-time migration [10], the source wavefield and receiver
wavefield are extrapolated forward and backward in time, respectively. Both
the source and receiver wavefields are decomposed into P- and S-wave
components, and various vector-based imaging conditions [33, 34, 35] can be
applied to generate images. We apply the 2D dot-product imaging condition

$I_{PP}=\sum_{t=0}^{T}\Big{(}U_{src\\{x\\}}^{P}U_{rec\\{x\\}}^{P}+U_{src\\{z\\}}^{P}U_{rec\\{z\\}}^{P}\Big{)}$
(10)

for the PP image and

$I_{PS}=\sum_{t=0}^{T}\Big{(}U_{src\\{x\\}}^{P}U_{rec\\{x\\}}^{S}+U_{src\\{z\\}}^{P}U_{rec\\{z\\}}^{S}\Big{)}$
(11)

for the PS image, where $T$ is the number of discrete time steps during
wavefield extrapolation; $U_{src}^{P}$, $U_{rec}^{P}$ and $U_{rec}^{S}$ are
the decomposed P-waves from the source wavefield, and decomposed P- and
S-waves from the receiver wavefield, respectively. All the wavefields in
equations 10 and 11 are understood to be functions of spatial position and
time $t$; these variables are omitted for simplicity in all the wavefield
expressions. Elastic RTMs are performed with DP and CNN-SF for PS wavefield
decomposition.

The Pluto model [36] is used in this test (Figure 14). The original Pluto
model has only the P-wave velocity defined. We approximate the S-wave velocity
model by multiplying the P-wave values by 0.6 at each grid point, and the
density has a constant value of 2.1 $\mathrm{g/cm^{3}}$.

Figure 14:  The Pluto P-wave velocity model. The S-wave velocity model is
approximated by multiplying the P-wave values by 0.6 at each grid point, and
the density has a constant value of 2.1 $\mathrm{g/cm^{3}}$.

The grid increments are h = 10 m in both spatial dimensions and time increment
dt = 1 ms. One hundred explosive sources with 15 Hz dominant frequency are
initiated sequentially from (x, z) = (0.5, 0.0) km to (13.5, 0.0) km with 0.13
km separation in the x direction. 1361 receivers are placed from (x, z) =
(0.0, 0.0) to (13.6, 0.0) km with 0.01 km separation. The stacked PP and PS
images using DP for PS decomposition are shown in Figure 15a and 15b,
respectively. The PS separation artifacts are visible along the high velocity
contrasts in the zoomed-in sections of Figure 15c and d. Compare with Figures
16a-16d, which are the PP, PS and zoomed-in images using CNN-SF for PS
decomposition; the image boundaries are more coherent than in Figure 15 as the
result of applying CNN-SF for PS decomposition.

Figure 15:  The (a) PP and (b) PS images of elastic RTM using DP for PS
decomposition; (c) and (d) are zoomed-in sections of boxed areas in (a) and
(b), respectively. Artifacts that are generated by DP remain in the stacked
images.  Figure 16:  The (a) PP and (b) PS images of elastic RTM using CNN-SF
for PS decomposition; (c) and (d) are zoomed-in sections of boxed areas in (a)
and (b), respectively. Compared with Figure 15, the artifacts are much
reduced.

##  VI Discussion

In the above tests, the dimensions of the CNN-tuned filters is fixed to be 15
$\times$ 15\. The accuracy of decomposition can be further improved by using
larger filters at the cost of computation efficiency. The ground truth
snapshots for training the CNN are the P-waves decomposed in the wavenumber
domain, thus the accuracy of CNN-tuned filters can never surpass that of the
wavenumber domain algorithm.

The computational complexity for PS decomposition using the 15 $\times$ 15
spatial filter set is $225N$, while the cost for the wavenumber domain
algorithm is $4N$log$(N)$, thus the latter is also much more efficient than
the former for small models. However, wavenumber domain algorithms require
Fourier transforms which consume large RAM and have difficulty being applied
locally in the model, while the spatial filters are HPC-friendly and can be
easily applied in parallel using multiple cores. Spatial filters are suitable
for processing target-oriented data sets, where PS decomposition can be
performed locally on a portion of the wavefield and thus further reducing the
cost.

The idea of using CNN-tuned spatial filters for PS decomposition can be
extended to anisotropic media. In that case, the spatial filters are non-
stationary across different positions within the model, because in anisotropic
wavefields, the polarization directions of P- and S-waves depend on the model
parameters [13], and several sets of filters can be tuned separately by
training with parameters from representative points in the model. The proposed
CNN-tuned spatial filters can also be extended to 3D, for which the number of
filters increases to 6 ($L_{x}$, $L_{y}$, $L_{z}$, $L_{xy}$, $L_{xz}$ and
$L_{yz}$) which need to be tuned as 3D convolutional operators.

##  VII Conclusions

We propose network-trained spatial filters for P- and S-wave decomposition in
isotropic elastic wavefields. The accuracy of decomposition using spatial
filters is improved after CNN-tuning with no extra computation costs.
Migration tests with synthetic data show that the proposed spatial filters do
not generate artifacts along high velocity contrasts as in the decoupled
propagation method. CNN-tuned spatial filters are suitable for parallel
computation and target-oriented imaging.

## Acknowledgment

The research leading to this paper is supported by the National Natural
Science Foundation of China (grant number: NSFC 41804108). The participation
of G.M. was supported by the UT-Dallas Geophysical Consortium.

## References

  * [1] X. Li, “Fracture detection using P-P and P-S waves in multicomponent sea-floor data,” in _Expanded Abstracts_. SEG, 1998, pp. 2056–2059.
  * [2] S. Knapp, N. Payne, and T. Johns, “Imaging through gas clouds: A case history from the Gulf of Mexico,” in _Expanded Abstracts_. SEG, 2001, pp. 3109–3113.
  * [3] B. Artman, I. Podladtchikov, and B. Witten, “Source location using time-reverse imaging,” _Geophysical Prospecting_ , vol. 58, pp. 861–873, 2010.
  * [4] J. Yang and H. Zhu, “Locating and monitoring microseismicity, hydraulic fracture and earthquake rupture using elastic time-reversal imaging,” _Geophysical Journal International_ , vol. 216, pp. 726–744, 2018.
  * [5] J. T. Kuo and F. T. Dai, “Kirchhoff elastic wave migration for the case of noncoincident source and receiver,” _Geophysics_ , vol. 49, pp. 1223–1238, 1984.
  * [6] F. T. Dai and J. T. Kuo, “Real data results of Kirchhoff elastic wave migration,” _Geophysics_ , vol. 51, pp. 1006–1011, 1986.
  * [7] K. Hokstad, “Multicomponent Kirchhoff migration,” _Geophysics_ , vol. 65, pp. 861–873, 2000.
  * [8] S. H. Gray, J. Etgen, J. Dellinger, and D. Whitmore, “Seismic migration problems and solutions,” _Geophysics_ , vol. 66, pp. 1622–1640, 2001.
  * [9] W. Chang and G. A. McMechan, “Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition,” _Geophysics_ , vol. 51, pp. 67–84, 1986.
  * [10] ——, “3-D elastic prestack reverse-time depth migration,” _Geophysics_ , vol. 59, pp. 597–609, 1994.
  * [11] N. D. Whitmore, “An imaging hierarchy for common-angle seismograms,” Ph.D. dissertation, The University of Tulsa, 1995.
  * [12] C. P. A. Wapenaar and G. C. Haimé, “Elastic extrapolation of primary seismic P- and S-waves,” _Geophysical Prospecting_ , vol. 38, pp. 23–60, 1990.
  * [13] J. Yan and P. Sava, “Elastic wave-mode separation for VTI media,” _Geophysics_ , vol. 74, pp. WB19–WB32, 2009.
  * [14] K. Aki and P. G. Richards, _Quantitative seismology, theory and methods_. W. H. Freeman and Co, 1980.
  * [15] R. Sun and G. A. McMechan, “Scalar reverse-time depth migration of prestack elastic seismic data,” _Geophysics_ , vol. 66, pp. 1518–1527, 2001.
  * [16] D. Ma and G. Zhu, “P- and S-wave separated elastic wave equation numerical modeling (in Chinese),” _Oil Geophysical Prospecting_ , vol. 38, pp. 482–486, 2003.
  * [17] Q. Zhang and G. A. McMechan, “2D and 3D elastic wavefield vector decomposition in the wavenumber domain for VTI media,” _Geophysics_ , vol. 75, pp. D13–D26, 2010.
  * [18] W. Wang, G. A. McMechan, and Q. Zhang, “Comparison of two algorithms for isotropic elastic P and S vector decomposition,” _Geophysics_ , vol. 80, pp. T147–T160, 2015.
  * [19] H. Zhu, “Elastic wavefield separation based on the Helmholtz decomposition,” _Geophysics_ , vol. 82, pp. S173–S183, 2017.
  * [20] X. Xiao and W. S. Leaney, “Local vertical seismic profiling (VSP) elastic reverse-time migration and migration resolution: Salt-flank imaging with transmitted P-to-S waves,” _Geophysics_ , vol. 75, pp. S35–S49, 2010.
  * [21] Y. Wei, H. Fu, Y. LI, and J. Yang, “A new P-wave reconstruction method for VSP data using conditional generative adversarial network,” in _Expanded Abstracts_. SEG, 2019, pp. 2528–2532.
  * [22] Y. Wei, Y. Li, J. Yang, J. Zong, J. Fang, and H. Fu, “Multi-task learning based P/S wave separation and reverse time migration for VSP,” in _Expanded Abstracts_. SEG, 2020, pp. 1671–1675.
  * [23] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel, “Backpropagation applied to handwritten zip code recognition,” _Neural Computation_ , vol. 1, pp. 541–551, 1989.
  * [24] H. Kaur, S. Fomel, and N. Pham, “Elastic wave-mode separation in heterogeneous anisotropic media using deep learning,” in _Expanded Abstracts_. SEG, 2019, pp. 2654–2658.
  * [25] J. Cheng and S. Fomel, “Fast algorithms for elastic-wave-mode separation and vector decomposition using low-rank approximation for anisotropic media,” _Geophysics_ , vol. 79, pp. C97–C110, 2014.
  * [26] W. Wang and G. A. McMechan, “Vector domain P and S decomposition in viscoelastic media,” in _Expanded Abstracts_. SEG, 2015, pp. 2153–2158.
  * [27] W. Wang, B. Hua, G. A. McMechan, and B. Duquet, “P and S decomposition in anisotropic media with localized low-rank approximations,” _Geophysics_ , vol. 83, pp. C13–C26, 2018.
  * [28] J. Paffenholz, B. McLain, J. Zaske, and P. J. Keliher, “Subsalt multiple attenuation and imaging: Observations from the Sigsbee2B synthetic dataset,” in _Expanded Abstracts_. SEG, 2002, pp. 2122–2125.
  * [29] J. Virieux, “P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method,” _Geophysics_ , vol. 51, pp. 889–901, 1986.
  * [30] D. Komatitsch and R. Martin, “An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation,” _Geophysics_ , vol. 72, pp. SM155–SM167, 2007.
  * [31] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” p. arXiv:1412.6980, 2014.
  * [32] G. S. Martin, K. J. Marfurt, and S. Larsen, “Marmousi-2: An updated model for the investigation of AVO in structurally complex areas,” in _Expanded Abstracts_. SEG, 2002, pp. 1979–1982.
  * [33] W. Wang and G. A. McMechan, “Vector-based elastic reverse-time migration,” _Geophysics_ , vol. 80, pp. S245–S258, 2015.
  * [34] Q. Du, G. Guo, Q. Zhao, X. Gong, C. Wang, and X. Li, “Vector-based elastic reverse time migration based on scalar imaging condition,” _Geophysics_ , vol. 82, pp. S111–S127, 2017.
  * [35] W. Wang and G. A. McMechan, “3D isotropic elastic reverse time migration using signed magnitudes of elastic data components,” _Journal of Seismic Exploration_ , vol. 28, 2019.
  * [36] D. Stoughton, J. Stefani, and S. Michell, “2D elastic model for wavefield investigations of subsalt objectives, deep water Gulf of Mexico.” EAGE, 2001, pp. A033–A036.