# Dephasing and Decorrelation of Spins in a Disordered Environment

Shreyas Raman  Department of Physics, Indian Institute of Science, Bangalore
560 012, India  Subroto Mukerjee  Department of Physics, Indian Institute of
Science, Bangalore 560 012, India

###### Abstract

Dephasing of spins is a major roadblock to scaling up the size of quantum
computing systems. We explore the possibility of utilizing highly disordered
environments which are in the Many-Body Localized phase to arrest this
dephasing. We embedded 2 ‘special’ spins in such a highly disordered
environment of Heisenberg spins to act as the target qubits and use the long-
time value of the spin-spin correlator
$\langle\vec{\sigma}_{i}\cdot\vec{\sigma}_{j}\rangle$ as an order parameter to
quantify the transition between the thermal and MBL phases of this system. It
is seen that the dephasing between spins, as encoded in this correlator, is
impeded in a disordered environment when the system is fully localized. The
order parameter yields a critical exponent, to characterize the transition
between the thermal and MBL phases, that appears to be robust to changes in
microscopic parameters of the system or the choice of pair of spins.

##  I Introduction

Quantum computation, which promises to be faster than its classical analog in
a range of applicationsGrover (1997) is based on the ability to store and
process information encoded in a collection of qubits, which are in
superpositions of states. An obstacle hindering the experimental
implementation of a complex quantum computational algorithm is quantum
decoherence, arising from the inevitable interaction between the system and
its environmentSackett  _et al._ (2000); Scully and Zubairy (1997). Further,
entanglement of qubits via quantum gatesSchmidt-Kaler  _et al._ (2003); Li
_et al._ (2003), which is an essential component of quantum computation,
results in faster loss of computational fidelityYu and Eberly (2004). Many
quantum error-correction protocols require the decoherence to be below a
certain thresholdZurek (2003). Thus, understanding and suppressing decoherence
in qubits is of paramount importance to the field of quantum computation.

The problems of single-qubit decoherence and its dynamical control have been
investigated thoroughlyKofman and Kurizki (2004); Greilich  _et al._ (2006)
and extended to multipartite decoherence controlGordon and Kurizki (2007). The
specific strategy to protect quantum information from decoherence depends on
the properties of the noise. When the coupling of different qubits is
uncorrelated, such as in the case of superconducting circuits, an effective
strategy is dynamical decoupling (DD), achieved by periodically applying a
series of pulses to the system to refocus the system-environment
evolutionViola and Lloyd (1998); Khodjasteh and Lidar (2005); Uhrig (2007).
Recent experiments have demonstrated significant improvement of coherence
times resulting from pulsed DDDamodarakurup  _et al._ (2009); de Lange  _et
al._ (2010); Ryan  _et al._ (2010); Souza  _et al._ (2011); Naydenov  _et al._
(2011). However, implementation of this method is challenging for systems with
fluctuating noise. Further, combining the decoupling scheme with gate
operations so that decoherence is suppressed is challenging.

The above mentioned methods have primarily focused on tuning the system under
consideration in an attempt to suppress decoherence. On the other hand, there
exists the possibility of modifying the environment parameters to do the same.
In this work, we propose tuning the environment to suppress decoherence by
introducing disorder into the environment. The motivation for this strategy
comes from noting that quantum systems reach thermal equilibrium starting from
an initial state through decoherenceHuse  _et al._ (2014). In such systems, we
typically find that the microscopic information associated with any initial
state is lost as the system evolves towards equilibrium. Information about any
initial state spreads over the entire system, rendering local measurements
incapable of any useful interpretationKim and Huse (2013); Deutsch (1991);
Srednicki (1994). Disorder-induced localization is a phenomenon that was first
shown by Anderson in 1958 in non-interacting systemsAnderson (1958). This
localized phase has been observed for systems of non-interacting phonons,
photons and matter wavesWeaver (1990); Schwartz  _et al._ (2007). Recent work
further shows that localization persists even in the presence of interactions
between particles, establishing the Many-Body-Localized (MBL) phase as a
robust, non-ergodic phase of quantum matter at finite temperaturesBasko  _et
al._ (2006); Gornyi  _et al._ (2005); Imbrie (2016a); Schreiber  _et al._
(2015a); Kondov  _et al._ (2015a); Fleishman and Anderson (1980); Oganesyan
and Huse (2007); Žnidarič  _et al._ (2008); Pal and Huse (2010); Bardarson
_et al._ (2012); Serbyn  _et al._ (2013a, b, 2014a, 2014b); Choi  _et al._
(2015); Yao  _et al._ (2015); Vosk and Altman (2014); Kjäll  _et al._ (2014);
Pekker  _et al._ (2014); Huse  _et al._ (2014); Luitz  _et al._ (2015);
Vasseur  _et al._ (2015); Vosk  _et al._ (2015); Agarwal  _et al._ (2015);
Devakul and Singh (2015); Schreiber  _et al._ (2015b); Imbrie (2016b); Kondov
_et al._ (2015b); yoon Choi  _et al._ (2016). Strong disorder, giving rise to
localization can prevent thermalization, leading to “memory” of initial states
even at late times. From an information theoretical perspective, the ability
of localized systems to avoid thermalization suggests their potential use as
memory resources. The presence of strong interactions leads to a logarithmic
growth of entanglement entropy, manifesting as an intrinsic, slow decoherence
mechanismBardarson  _et al._ (2012); Serbyn  _et al._ (2013a, b, 2014a,
2014b); Huse  _et al._ (2014), which has been experimentally observedLukin
_et al._ (2019); Brydges  _et al._ (2019). A recent numerical study of the
dynamics on the ergodic side of the transition in a Floquet model of many-body
localization in large systems sizes has revealed clear slow dynamics for a
given system size. This is reflected in a stretched exponential decay of the
autocorrelation function(equal position, unequal time correlator) and
subballistic spreading of the entanglement entropy; with the dynamics
consistently speeding up with increasing system size and decreasing disorder
strengthLezama  _et al._ (2019).

The ergodic to MBL transition has been a topic of intense study for many
years. However, the exact nature of the transition is still not well
understoodAlet and Laflorencie (2018); Khemani  _et al._ (2017); Morningstar
_et al._ (2022); Crowley and Chandran (2020). A part of the problem lies in
the fact that it appears not to be describable in terms of symmetry breaking
with an appropriate order parameter like a conventional continuous transition.
Thus, it is not obvious whether the transition can be characterized by
critical exponents (and if so, by how many) that are robust to microscopic
details. In fact, the ergodic to MBL transition is “dynamic” in the sense that
it separates a phase (the ergodic phase) which dynamically equilibrates
starting from an arbitrary initial condition from one which does not (the MBL
phase). Since the dynamics is completely determined by the spectrum of energy
eigenvalues and eigenstates, the ergodic to MBL transition is captured in the
change of properties of the eigenstates and the distribution of eigenvalues
across the transition. Even in the absence of any obvious symmetry breaking,
attempts have been made to characterize the transition in terms of different
order parameters. Examples of such order parameters are the scaling of the
entanglement entropy with system sizeAlet and Laflorencie (2018), the
eigenstate to eigenstate fluctuation of physical observablesDeng  _et al._
(2017); Li  _et al._ (2015), the density imbalance at long times starting from
an initial non-zero valueyoon Choi  _et al._ (2016), the fluctuation of
entanglement entropy from eigenstate to eigenstate, one disorder realization
to another and from cut to cut to within a single systemKhemani  _et al._
(2017). However, the nature of the transition and even whether it is
continuous or not remains inconclusive from these studies.

Several studies have also focused on obtaining a diverging correlation length
to characterize the transition. Studies on finite-sized systemsKjäll  _et al._
(2014); Luitz  _et al._ (2015); Khemani  _et al._ (2017) of microscopic models
with quenched randomness have produced correlation length exponents $\nu$ that
violate the Chayes-Chayes-Fisher-Spencer (CCFS) bound $\nu>2$Li  _et al._
(2015) suggesting that they are not in the proper scaling regime to obtain an
accurate value of $\nu$. A complementary approach has utilized renormalization
group methods for phenomenological models with coarse-grained interacting
thermal and insulating regions Vosk  _et al._ (2015); Goremykina  _et al._
(2019); Thiery  _et al._ (2018). The most fruitful of these has been based on
the avalanche mechanism Thiery  _et al._ (2018) have produced exponents
consistent with the CCFS bound and in fact, even $\nu=\infty$ suggesting that
the transition is Kosterlitz-Thouless likeDumitrescu  _et al._ (2019).

Numerics on finite-sized systems have so far not found rare thermal regions
required to seed avalanches. This has been attributed to the insufficient
sizes of these systems. It has been argued that the “true” MBL to thermal
transition in the thermodynamic limit which occurs via avalanches is at a
disorder strength significantly larger than obtained from finite-size
numericsMorningstar  _et al._ (2022). Additionally, the “transition” observed
in numerical studies has been claimed to be a crossover (albeit a fairly sharp
one) from a finite-size MBL phase destabilized by system-wide many-body
resonances to a thermal phaseMorningstar  _et al._ (2022); Crowley and
Chandran (2020). It should be emphasized that the true MBL transition remains
out of reach numerically and also experimentally owing to limitations of
system size and time scales. Thus, it is only the finite-size crossover, which
can be studied in microscopic models without an appeal to phenomenology and it
is still not well understood. The existence of the MBL phase in the
thermodynamic limit has also been questioned with no proper evidence currently
available to resolve the issue Šuntajs  _et al._ (2020a, b); Sels and
Polkovnikov (2021).

In this paper, we focus on the sharp MBL to thermal crossover accessible to
numerics on finite-sized systems. To be consistent with the prevailing
terminology, we refer to this crossover as a “transition” with the
understanding that it is not expected to be one in the thermodynamic limit.
The particular quantity we study is the equal-time correlator of a pair of
“special” spins embedded in the rest of of the system, which can be in the
thermal or the MBL phase depending on the strength of the disorder. These
special spins, which are not subjected to disorder themselves, interact with
the rest of the system. Our system is thus a very simple model of two qubits
in an environment whose properties can be tuned with changing disorder
strength. We show that the asymptotic value of the above-mentioned correlator
is suppressed with increasing disorder strength. Further, we argue that this
quantity can act as an “order parameter” for the transition and numerically
extract a critical exponent that appears to be fairly robust to changes in the
microscopic parameters.

This paper is organized as follows : In section II, we introduce the system
under consideration along with its ratio of successive level spacingsAlet and
Laflorencie (2018); Oganesyan and Huse (2007), to identify the transition
between the thermal and MBL phase. Further, we define the spin-spin
correlator, which we employ as a diagnostic to quantify decoherence. In
section III, we demonstrate that the correlator serves as an order parameter
for the transition and obtain a critical exponent from the same for the
transition for different system parameters. The data obtained has been
tabulated in Section IV.

##  II Spin-Spin Correlator in a Disordered Environment

Consider the following one-dimensional lattice model, the spin-$\frac{1}{2}$
Heisenberg chain in a random magnetic field

$\hat{\mathcal{H}}=\sum_{i}\left(\hat{\sigma}_{i}^{x}\hat{\sigma}_{i+1}^{x}+\hat{\sigma}_{i}^{y}\hat{\sigma}_{i+1}^{y}+\hat{\sigma}_{i}^{z}\hat{\sigma}_{i+1}^{z}\right)+\sum_{i}h_{i}\hat{\sigma}_{i}^{z}$
(1)

where the system has periodic boundary conditions, $\sigma^{x,y,z}$ are the
Pauli spin matrices and the $h_{i}$ are drawn from a uniform distribution
$\left[-h,h\right]$. We first use the Jordan-Wigner transformation to convert
this to a model of hard-core bosons for computational fidelity. With some re-
scaling, this can be written as a simple tight-binding system of hard-core
bosons with nearest neighbour hopping and interaction with random on-site
potential, i.e

$\hat{\mathcal{H}}=\sum_{i}\left(\hat{b}_{i}^{\dagger}\hat{b}_{i+1}+\hat{b}_{i+1}^{\dagger}\hat{b}_{i}+2\hat{n}_{i}\hat{n}_{i+1}\right)+\sum_{i}h_{i}\hat{n}_{i}$
(2)

where $b$ and $b^{\dagger}$ are the hard-core boson annihilation and creation
operators respectively and $n_{i}$ is the occupation of site $i$. Note that
this Hamiltonian preserves boson number, just like the original Hamiltonian
preserves total spin. The numerical analysis of this paper will be done on a
system with equal number of up and down spins, or in terms of bosons, a system
with half-filling.

Eigenvalue Statistics : In the thermal phase, the eigenvalues follow the
statistics of Random Matrix Theory. As the Hamiltonian has only real entries,
this will be the Gaussian Orthogonal Ensemble (GOE). A useful computational
quantifier to verify this is the ratio of consecutive level spacings Alet and
Laflorencie (2018); Oganesyan and Huse (2007)

$r_{n}=\frac{\min\left(E_{n}-E_{n-1},E_{n-1}-E_{n-2}\right)}{\max\left(E_{n}-E_{n-1},E_{n-1}-E_{n-2}\right)}$
(3)

Its average value in the GOE can be computed as $\langle
r\rangle_{\text{GOE}}\approx 0.5307$. On the other hand, in the MBL phase, the
eigenvalues are non-correlated and hence they follow Poisson statistics. The
average value in the Poisson distribution can be computed as $\langle
r\rangle_{\text{Poisson}}\approx 0.386$. For this system, we compute the
average ratio by averaging over disorder realizations.

###  II.1 Modified System with special spins

We explore the viability of the spin-spin correlator between 2 ‘special’ spins
as an order parameter for the phase transition between the thermal and MBL
phases. The 2 special spins are set such that they do not have a random on-
site potential, while every other spin does. In terms of bosons, this is
equivalent to picking 2 special sites($i$ and $j$), which do not have any
random on-site potential. The Hamiltonian is :

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=\sum_{k}\left(\hat{b}_{k}^{\dagger}\hat{b}_{k+1}+\hat{b}_{k+1}^{\dagger}\hat{b}_{k}+2\hat{n}_{k}\hat{n}_{k+1}\right)+\sum_{k;k\neq
i,j}h_{k}\hat{n}_{k}$ (4) $\displaystyle\hat{\mathcal{H}}(V)$
$\displaystyle=\sum_{k}\left(\hat{b}_{k}^{\dagger}\hat{b}_{k+1}+\hat{b}_{k+1}^{\dagger}\hat{b}_{k}+V\hat{n}_{k}\hat{n}_{k+1}\right)+\sum_{k;k\neq
i,j}h_{k}\hat{n}_{k}$ (5)

where $V$ is a parameter introduced as a microscopic system parameter - the
Interaction Potential Strength - that can be tuned (explained later). When not
mentioned, $V=2$.

Comparing the ratio of successive levels of the original system with the
modified system (FIG. 1) confirms that in terms of level-spacing statistics,
there is no major change in the system. The spin-spin correlator we want to
evaluate is $\langle\vec{\sigma_{i}}\cdot\vec{\sigma}_{j}\rangle$. In terms of
boson operators, this is :

$\langle\vec{\sigma_{i}}\cdot\vec{\sigma_{j}}\rangle\equiv
4n_{i}n_{j}-2n_{i}-2n_{j}+1+2\langle b_{i}^{\dagger}b_{j}\rangle+2\langle
b_{j}^{\dagger}b_{i}\rangle$ (6)

Figure 1: Ratio of successive levels averaged over 2000 disorder realizations.
Ratio plotted for different choices of special spins (including no special
pair). Inset : Ratio of successive levels for a 14 site system for the spin
pair 3, 8. Figure 2: Correlator vs. Time for $d=5$ and for different values of
disorder amplitude; averaged over 2000 disorder realizations and 210 initial
configurations. Inset : Correlator vs Time for 14 Site system for $h=2,6,10$,
averaged over 1500 disorder realizations and 792 initial configurations

The system is initialized such that at time $t=0$, it is in a product state
where both the special spins are initialized to ‘up’. We then average over all
possible initial product states that have this property. This ensures the
initial correlator is exactly 1. The correlator as a function of time was
plotted for the system, averaged over disorder realizations and initial
configurations(FIG. 2). A fixed distance between spin pairs was chosen,
defined as $d\equiv|j-i|$.

The trend shown by the long-time value of the correlator indicates that it
could be a good order parameter for characterizing the phase transition
between the thermal and MBL phases. For small disorder strength $h$, the long-
time correlator tends to the same value, characteristic of the system
thermalizing. In the infinite system size limit, this value tends to 0. From
here, as $h$ increases, the long-time correlator also increases. The trend for
the long-time correlator as a function of disorder strength is shown in FIG.
3.

Figure 3: Long-Time value of the Correlator vs. Disorder Amplitude for a 12
site system; for different values of interaction potential and choice of spin
pair. The trend shown by the long-time correlator motivates its use as an
order parameter to quantify this transition and subsequent extraction of a
power law close to the transition point.

##  III Transition Points and Critical Exponents

To find the transition point of an $N$ site system, we found the intersection
point of the ratio of successive levels between an $N$ site system and an
$N-2$ site system(FIG. 4). All the other parameters of the system were kept
the same.

A power law was extracted from the trend of long-time correlator vs. disorder
strength close to the transition point. To do this, we assume the relation is
of the form (to the right of the transition point) :

$\langle\vec{\sigma}_{i}-\vec{\sigma}_{j}\rangle-\langle\vec{\sigma}_{i}-\vec{\sigma}_{j}\rangle_{c}=A(h-h_{c})^{\beta}$
(7)

where the subscript $c$ stands for the value at the critical (transition)
point.

The robustness of this power-law exponent can be checked against two major
parameters :

  1. 1.

Choice of spin pair

  2. 2.

Interaction Potential

###  III.1 Choice of Spin Pair

As shown by FIG 1, the physics of the system does not depend on the choice of
spin pair and transition point of the system is independent of choice of spin
pair.

###  III.2 Interaction Potential

Figure 4: Ratio of successive levels showing the transition points for
different values of interaction potential. Inset : Ratio of successive levels
showing the transition points for a 14 site system. Figure 5: Close to the
transition point, it was assumed that the long-time correlator follows a power
law as a function of disorder strength of the form
$\langle\vec{\sigma}_{i}\cdot\vec{\sigma}_{j}\rangle-\langle\vec{\sigma_{i}}\cdot\vec{\sigma}_{j}\rangle_{c}\sim(h-h_{c})^{\beta}$.
The plot shows the fitted Critical Exponent close to the transition point for
different values of interaction potential and choice of spin pair. The
exponent obtained has an average value that seems to tend towards $1$ with
increasing system size and is robust to system parameters (choice of spin pair
and interacting potential).

This far in the paper, we have mainly discussed $V=2$. As we explore the
robustness of the long-time value of the correlator as an order parameter, we
turn to changing $V$ as a microscopic parameter that changes the physics of
the system. As the form of the Hamiltonian changes on changing the interaction
potential $V$, the transition point of the system also changes. In order to
obtain the new transition point, we plot the ratio of successive levels for 2
system sizes with the changed potential and locate the new intersection point
(FIG. 4).

A power law is fitted to the long-time value of the correlator to the right of
the transition points as obtained from the intersection in FIG. 4. The
exponent values for some values of $V$ are obtained in FIG. 5 and the data has
been tabulated in Table 1 and Table 2.

##  IV Data

Table 1: Exponent $\beta$ cross-referenced by Potential and distance between spin pairs for 12 Site system $V$ $d$ | 4 | 5 | 6 | Average
---|---|---|---|---
0.5 | 1.1 | 1.15 | 1.15 |  1.13 $\pm$0.03
1.0 | 1.09 | 1.08 | 1.09 |  1.09 $\pm$0.01
1.5 | 1.07 | 1.12 | 1.04 |  1.08 $\pm$0.04
2.0 | 1.11 | 1.12 | 1.05 |  1.09 $\pm$0.04
2.5 | 1.1 | 1.07 | 1.09 |  1.09 $\pm$0.01
3.0 | 1.02 | 1.09 | 1.18 |  1.1 $\pm$0.08
Average |  1.08 $\pm$ 0.03 |  1.11 $\pm$ 0.03 |  1.10 $\pm$0.05 |  1.10 $\pm$ 0.04
Table 2: Exponent $\beta$ for different values of interacting potential for 14 site system with $d=5$ Potential($V$) |  Exponent($\beta$)
---|---
0.5 | 1.03
1.0 | 1.00
1.5 | 1.05
2.0 | 1.05
2.5 | 1.00
3.0 | 0.99
Average | 1.02 $\pm$ 0.03

##  V Conclusion and Discussion

From the calculations performed on a modified Heisenberg spin chain with a
pair of special spins, it appears that the equal time spin-spin correlator of
these spins can serve as an order-parameter for the transition of the rest of
the system, from the thermal to MBL phase. The asymptotic value of this
quantity is non-zero in the MBL phase of the environment and equal to zero
(modulo finite-size effects) in the MBL phase. This clearly demonstrates that
the decoherence of these spins is impeded in the presence of a disordered
environment, which can find use in the design of quantum memories and quantum
computing platforms. Interestingly, we find that the correlator yields a
critical exponent that appears to have a value close to unity irrespective of
the choice of parameters of the Hamiltonian or the distance between the
special spins. This appears to justify the use of the correlator as an order
parameter and hints at underlying universal features of the transition. The
existence of such a universality would be interesting since, as mentioned in
the introduction, it has been recently argued that the transition being
studied here is, in reality only a sharp crossover. Whether the apparent
robustness of the exponent obtained survives to higher system sizes is
currently a matter of conjecture because of numerical limitations.
Nevertheless, our observation is still significant given that all existing
numerical studies on microscopic systems and indeed experimentally accessible
systems appear to only be able to probe the finite-sized MBL transition and
not the one that has been claimed to occur in the thermodynamic limit.

##  VI Acknowledgments

SM thanks the Department of Science and Technology, Govt. of India for funding
for this work through the Quantum Information Science and Technology (QuST)
scheme.

## References

  * Grover (1997) L. K. Grover, Phys. Rev. Lett. 79, 325 (1997).
  * Sackett  _et al._ (2000) C. A. Sackett, D. Kielpinski, B. E. King, C. Langer, V. Meyer, C. J. Myatt, M. Rowe, Q. A. Turchette, W. M. Itano, D. J. Wineland, and C. Monroe, Nature 404, 256 (2000).
  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Schmidt-Kaler  _et al._ (2003) Schmidt-Kaler, F. Häffner, H. Riebe, M. Gulde, S. Lancaster, G. P. T. Deuschle, T. Becher, C. Roos, C. F. Eschner, J. Blatt, and Rainer, Nature 422, 408 (2003).
  * Li  _et al._ (2003) X. Li, Y. Wu, D. Steel, D. Gammon, T. H. Stievater, D. S. Katzer, D. Park, C. Piermarocchi, and L. J. Sham, Science 301, 809 (2003), https://www.science.org/doi/pdf/10.1126/science.1083800 .
  * Yu and Eberly (2004) T. Yu and J. H. Eberly, Phys. Rev. Lett. 93, 140404 (2004).
  * Zurek (2003) W. H. Zurek, Rev. Mod. Phys. 75, 715 (2003).
  * Kofman and Kurizki (2004) A. G. Kofman and G. Kurizki, Phys. Rev. Lett. 93, 130406 (2004).
  * Greilich  _et al._ (2006) A. Greilich, D. R. Yakovlev, A. Shabaev, A. L. Efros, I. A. Yugova, R. Oulton, V. Stavarache, D. Reuter, A. Wieck, and M. Bayer, Science 313, 341 (2006), https://www.science.org/doi/pdf/10.1126/science.1128215 .
  * Gordon and Kurizki (2007) G. Gordon and G. Kurizki, Physical Review A 76 (2007), 10.1103/PhysRevA.76.042310.
  * Viola and Lloyd (1998) L. Viola and S. Lloyd, Phys. Rev. A 58, 2733 (1998).
  * Khodjasteh and Lidar (2005) K. Khodjasteh and D. A. Lidar, Phys. Rev. Lett. 95, 180501 (2005).
  * Uhrig (2007) G. S. Uhrig, Phys. Rev. Lett. 98, 100504 (2007).
  * Damodarakurup  _et al._ (2009) S. Damodarakurup, M. Lucamarini, G. Di Giuseppe, D. Vitali, and P. Tombesi, Phys. Rev. Lett. 103, 040502 (2009).
  * de Lange  _et al._ (2010) G. de Lange, Z. H. Wang, D. Ristè, V. V. Dobrovitski, and R. Hanson, Science 330, 60 (2010), https://www.science.org/doi/pdf/10.1126/science.1192739 .
  * Ryan  _et al._ (2010) C. A. Ryan, J. S. Hodges, and D. G. Cory, Phys. Rev. Lett. 105, 200402 (2010).
  * Souza  _et al._ (2011) A. M. Souza, G. A. Álvarez, and D. Suter, Phys. Rev. Lett. 106, 240501 (2011).
  * Naydenov  _et al._ (2011) B. Naydenov, F. Dolde, L. T. Hall, C. Shin, H. Fedder, L. C. L. Hollenberg, F. Jelezko, and J. Wrachtrup, Phys. Rev. B 83, 081201 (2011).
  * Huse  _et al._ (2014) D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014).
  * Kim and Huse (2013) H. Kim and D. A. Huse, Phys. Rev. Lett. 111, 127205 (2013).
  * Deutsch (1991) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Weaver (1990) R. Weaver, Wave Motion 12, 129 (1990).
  * Schwartz  _et al._ (2007) Schwartz, T. Bartal, G. Fishman, S. Segev, and Mordechai, Nature 446, 52 (2007).
  * Basko  _et al._ (2006) D. Basko, I. Aleiner, and B. Altshuler, Annals of Physics 321, 1126 (2006).
  * Gornyi  _et al._ (2005) I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. 95, 206603 (2005).
  * Imbrie (2016a) J. Z. Imbrie, Phys. Rev. Lett. 117, 027201 (2016a).
  * Schreiber  _et al._ (2015a) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842 (2015a), https://www.science.org/doi/pdf/10.1126/science.aaa7432 .
  * Kondov  _et al._ (2015a) S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco, Phys. Rev. Lett. 114, 083002 (2015a).
  * Fleishman and Anderson (1980) L. Fleishman and P. W. Anderson, Phys. Rev. B 21, 2366 (1980).
  * Oganesyan and Huse (2007) V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).
  * Žnidarič  _et al._ (2008) M. Žnidarič, T. c. v. Prosen, and P. Prelovšek, Phys. Rev. B 77, 064426 (2008).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Bardarson  _et al._ (2012) J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
  * Serbyn  _et al._ (2013a) M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 110, 260601 (2013a).
  * Serbyn  _et al._ (2013b) M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 111, 127201 (2013b).
  * Serbyn  _et al._ (2014a) M. Serbyn, M. Knap, S. Gopalakrishnan, Z. Papić, N. Y. Yao, C. R. Laumann, D. A. Abanin, M. D. Lukin, and E. A. Demler, Phys. Rev. Lett. 113, 147204 (2014a).
  * Serbyn  _et al._ (2014b) M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. B 90, 174302 (2014b).
  * Choi  _et al._ (2015) S. Choi, N. Y. Yao, S. Gopalakrishnan, and M. D. Lukin, “Quantum control of many-body localized states,”  (2015).
  * Yao  _et al._ (2015) N. Y. Yao, C. R. Laumann, and A. Vishwanath, “Many-body localization protected quantum state transfer,”  (2015).
  * Vosk and Altman (2014) R. Vosk and E. Altman, Phys. Rev. Lett. 112, 217204 (2014).
  * Kjäll  _et al._ (2014) J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 113, 107204 (2014).
  * Pekker  _et al._ (2014) D. Pekker, G. Refael, E. Altman, E. Demler, and V. Oganesyan, Phys. Rev. X 4, 011052 (2014).
  * Luitz  _et al._ (2015) D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 91, 081103 (2015).
  * Vasseur  _et al._ (2015) R. Vasseur, S. A. Parameswaran, and J. E. Moore, Phys. Rev. B 91, 140202 (2015).
  * Vosk  _et al._ (2015) R. Vosk, D. A. Huse, and E. Altman, Phys. Rev. X 5, 031032 (2015).
  * Agarwal  _et al._ (2015) K. Agarwal, S. Gopalakrishnan, M. Knap, M. Müller, and E. Demler, Phys. Rev. Lett. 114, 160401 (2015).
  * Devakul and Singh (2015) T. Devakul and R. R. P. Singh, Phys. Rev. Lett. 115, 187201 (2015).
  * Schreiber  _et al._ (2015b) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842 (2015b), https://www.science.org/doi/pdf/10.1126/science.aaa7432 .
  * Imbrie (2016b) J. Z. Imbrie, Phys. Rev. Lett. 117, 027201 (2016b).
  * Kondov  _et al._ (2015b) S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco, Phys. Rev. Lett. 114, 083002 (2015b).
  * yoon Choi  _et al._ (2016) J. yoon Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Science 352, 1547 (2016), https://www.science.org/doi/pdf/10.1126/science.aaf8834 .
  * Lukin  _et al._ (2019) A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, and M. Greiner, Science 364, 256 (2019), https://www.science.org/doi/pdf/10.1126/science.aau0818 .
  * Brydges  _et al._ (2019) T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Science 364, 260 (2019), https://www.science.org/doi/pdf/10.1126/science.aau4963 .
  * Lezama  _et al._ (2019) T. L. M. Lezama, S. Bera, and J. H. Bardarson, Physical Review B 99 (2019), 10.1103/physrevb.99.161106.
  * Alet and Laflorencie (2018) F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498 (2018).
  * Khemani  _et al._ (2017) V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse, Phys. Rev. X 7, 021013 (2017).
  * Morningstar  _et al._ (2022) A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Physical Review B 105 (2022), 10.1103/physrevb.105.174205.
  * Crowley and Chandran (2020) P. J. D. Crowley and A. Chandran, “A constructive theory of the numerically accessible many-body localized to thermal crossover,”  (2020).
  * Deng  _et al._ (2017) D.-L. Deng, S. Ganeshan, X. Li, R. Modak, S. Mukerjee, and J. H. Pixley, Annalen der Physik 529, 1600399 (2017), https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.201600399 .
  * Li  _et al._ (2015) X. Li, S. Ganeshan, J. H. Pixley, and S. Das Sarma, Phys. Rev. Lett. 115, 186601 (2015).
  * Goremykina  _et al._ (2019) A. Goremykina, R. Vasseur, and M. Serbyn, Phys. Rev. Lett. 122, 040601 (2019).
  * Thiery  _et al._ (2018) T. Thiery, F. m. c. Huveneers, M. Müller, and W. De Roeck, Phys. Rev. Lett. 121, 140601 (2018).
  * Dumitrescu  _et al._ (2019) P. T. Dumitrescu, A. Goremykina, S. A. Parameswaran, M. Serbyn, and R. Vasseur, Phys. Rev. B 99, 094205 (2019).
  * Šuntajs  _et al._ (2020a) J. Šuntajs, J. Bonča, T. c. v. Prosen, and L. Vidmar, Phys. Rev. E 102, 062144 (2020a).
  * Šuntajs  _et al._ (2020b) J. Šuntajs, J. Bonča, T. c. v. Prosen, and L. Vidmar, Phys. Rev. B 102, 064207 (2020b).
  * Sels and Polkovnikov (2021) D. Sels and A. Polkovnikov, Phys. Rev. E 104, 054105 (2021).


# Derivation of Analytic Formulas for the Sample Moments of the Sample
Correlation over Permutations of Data

Marc Jaffrey111University of Washington<EMAIL_ADDRESS>222RootFault LLC
Michael Dushkoff333Rochester Institute of Technology<EMAIL_ADDRESS>††footnotemark:

###### Abstract

Pearson’s correlation $\rho$ is among the mostly widely reported measures of
association. The strength of the statistical evidence for linear association
is determined by the p-value of a hypothesis test. If the true distribution of
a dataset is bivariate normal, then under specific data transformations a
t-statistic returns the exact p-value, otherwise it is an approximation.
Alternatively, the p-value can be estimated by analyzing the distribution of
the sample correlation, $r$, under permutations of the data. Moment
approximations of this distribution are not as widely used since estimation of
the moments themselves are numerically intensive with greater uncertainties.
In this paper we derive an inductive formula allowing for the analytic
expression of the sample moments of the sample correlation under permutations
of the data in terms of the central moments of the data. These formulas placed
in a proper statistical framework could open up the possibility of new
estimation methods for computing the p-value.

##  1 Introduction

Pearson’s product-moment correlation [1],

$\rho=\frac{COV(X,Y)}{\sigma_{X}\sigma_{Y}}$ (1)

the generalization of Galton’s regression coefficient, [2, 3] is perhaps the
most important measure of data association, orienting scientists in the
direction of discovery since its introduction. As an inferred measure of
linear association between two variables, it is a primary tool in statistical
analysis and the search for meaningful variable relationships, (see Figure 1).
The subject of its interpretation and nuances surrounding its use have been
extensively covered in the literature with still plenty of room for debate
today, [4, 5].

Given a dataset with two variables, in hypothesis testing,

$\displaystyle H_{0}:\rho=0$ (2) $\displaystyle H_{1}:\rho\neq 0$

the p-value of the test determines the strength of the statistical evidence
for whether to accept or reject $H_{0}$ [6]. Widely contentious, the use of
p-value for determining statistical significance is under critical debate, a
subject that is extensively addressed in the recent literature [7, 8, 9, 10,
11, 12]. Notwithstanding this controversy, p-value is still in widespread use
today, and as such, methods of its computation are relevant to discuss.

Presented is a branching inductive formula for computing the sample moments of
the sample correlation coefficient $r$,

$r=\frac{\sum_{i}(x_{i}-\mu_{x})(y_{i}-\mu_{y})}{\sqrt{\sum_{k}(x_{k}-\mu_{x})^{2}\sum_{j}(y_{j}-\mu_{y})^{2}}},$
(3)

over permutations of the data in terms of the central moments of the data,
denoted

$\langle r_{\Pi}^{k}\rangle=\frac{1}{|\Pi|}\sum_{\pi\in\Pi}r_{\pi}^{k}$ (4)

for $k\in\mathbb{N}$, where $\Pi$ is the set of all possible permutations on
the data and $r_{\pi}$ denotes the sample correlation over the permutation
$\pi\in\Pi$ of the data. These formulas are interesting by themselves and
placed in a proper statistical framework could potentially open the door to
the possibility of computationally efficient methods for computing the p-value
of a hypothesis test of $\rho$.

Figure 1:  All three synthetic examples clearly exhibit strong functional
relationships. Under a hypothesis test for the presence of a linear
relationship under $\rho$, the first two examples strongly reject the null
hypothesis, both with p-values $<\mathit{1\times 10^{-4}}$, while in the last
example the null hypothesis is accepted. Even though the middle example
reports a moderate linear association, it is clearly non-linear being a cubic
polynomial. Pearson’s correlation is perhaps better interpreted as a measure
of the suitability of a linear approximation to the data.

##  2 Main Result: Inductive Formula for $\langle r_{\Pi}^{k}\rangle$

###  2.1 Notation

Given a dataset $D_{n}=\\{(x_{i},y_{i})|\ i=1,..,n\\}\subset\mathbb{R}^{2}$,
let $\Pi=Perm(n)$ be the set of all permutations on $\\{1,\dots,n\\}$. For
$\pi=(\pi_{1},\dots,\pi_{n})\in\Pi$, defined such that
$i\underset{\pi}{\to}\pi_{i}$, define

$\pi(D_{n})=\\{(x_{i},y_{\pi_{i}})|i=1,..,n\\},$ (5)

where the permutation $\pi$ acts on the $y$-coordinate alone. Denote Pearson’s
sample correlation over a permutation $\pi$ of the data as
$r_{\pi}=r(\pi(D_{n}))$ and define

$\displaystyle r_{\Pi}:[-1,1]\to[0,1],$ (6)

the distribution of $r$ over the permutations of $D_{n}$. Lastly, let
$\hat{x}_{i}\equiv(x_{i}-\mu_{x})$ and $\hat{y}_{i}\equiv(y_{i}-\mu_{y})$.

###  2.2 Main Result

Given $D_{n}$, fix $k$. The $k^{th}$ moment of $r_{\Pi},$ is given by

$\langle
r_{\Pi}^{k}\rangle=\frac{1}{n^{k}n!\hat{\sigma}_{x}^{k}\hat{\sigma}_{y}^{k}}\sum_{m=1}^{k}\
\
\sum_{n_{1}+...+n_{m}=k}\binom{k}{n_{1},..,n_{m}}^{*}X_{(n_{1},..,n_{m})}^{m,k}(n-m)!Y_{(n_{1},..,n_{m})}^{m,k}\cdot
h_{n,m}$ (7)

$X_{(n_{1},..,n_{m})}^{m,k}\text{ and }Y_{(n_{1},..,n_{m})}^{m,k}$ are
branching inductive formulas, such that for $1<m\leq k$ and for fixed non-zero
positive integers $n_{1}+...+n_{m}=k$ we have:

$X_{(n_{1},..,n_{m})}^{m,k}=\ n\langle\hat{x}^{n_{m}}\rangle
X_{(n_{1},...,n_{m-1})}^{m-1,k-n_{m}}-\sum_{j}X_{(n_{1}+n_{m}\delta_{ij},...,n_{m-1}+n_{m}\delta_{(m-1)j})}^{m-1,k}$
(8)

where $\delta_{ij}$ is the standard delta function, and similarly defined

$Y_{(n_{1},..,n_{m})}^{m,k}=\ n\langle\hat{y}^{n_{m}}\rangle
Y_{(n_{1},...,n_{m-1})}^{m-1,k-n_{m}}-\sum_{j}Y_{(n_{1}+n_{m}\delta_{ij},...,n_{m-1}+n_{m}\delta_{(m-1)j})}^{m-1,k}$
(9)

For $m=1$, we have

$X_{(k)}^{1,k}=n\langle\hat{x}^{k}\rangle\text{ and
}Y_{(k)}^{1,k}=n\langle\hat{y}^{k}\rangle$ (10)

The stared multinomial coefficient

$\binom{k}{n_{1},\dots,n_{m}}^{*}=\frac{1}{d_{1}!\cdots
d_{r}!}\binom{k}{n_{1},\dots,n_{m}}$ (11)

is an adjustment of the usual multinomial coefficient accounting for
degeneracy in $n_{1},...,n_{m}$. Partitioning $n_{1},...m_{m}$ into subsets,
$g_{1},..,g_{r}$, by the equivalence relation $n_{i}\equiv n_{j}\iff
n_{i}=n_{j}$, then $d_{i}=|g_{i}|$. This degeneracy leads to an over counting
represented by the multinomial coefficient which is correct by dividing out
$d_{1}!\cdots d_{r}!$.

The term $h_{n,m}$ accounts for the inability to compute higher order terms in
the sum when the number of data points is less than the moment order being
computed, by setting them to zero in the formula,

$h_{n,m}=\left\\{\begin{array}[]{rcl}0&\mbox{for}&n-m<0\\\
1&\mbox{for}&n-m\geq 0\end{array}\right.$ (12)

Lastly, for notational simplicity, in (7) the following convention is
employed,

$\hat{\sigma}_{z}=\sqrt{\frac{1}{n}\sum(z_{i}-\mu_{z})^{2}}$ (13)

##  3 Derivation

Starting from (4) and (1), we have

$\displaystyle\langle r_{\Pi}^{k}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{n!}\sum_{\pi\in\Pi}r_{\pi}^{k}$ (14) $\displaystyle=$
$\displaystyle\frac{1}{n!}\sum_{\pi\in\Pi}\Bigg{[}\frac{\sum_{i}(x_{i}-\mu_{x})(y_{\pi_{i}}-\mu_{y})}{\sqrt{\sum_{k}(x_{k}-\mu_{x})^{2}\sum_{j}(y_{\pi_{j}}-\mu_{y})^{2}}}\Bigg{]}^{k}$
(15) $\displaystyle=$
$\displaystyle\frac{1}{n^{k}n!\hat{\sigma}_{x}^{k}\hat{\sigma}_{y}^{k}}\sum_{\pi\in\Pi}\sum_{i_{i},..,i_{k}}\hat{x}_{i_{1}}...\hat{x}_{i_{k}}\hat{y}_{\pi_{i_{1}}}...\hat{y}_{\pi_{i_{k}}}$
(16)

Reorganizing the summation over the indices, let $1\leq m\leq k$ be the number
of distinct indices and the non-zero positive integers $n_{1}+\cdots+n_{m}=k$
their multiplicity, then (16) becomes

$\displaystyle\frac{1}{n^{k}n!\hat{\sigma}_{x}^{k}\hat{\sigma}_{y}^{k}}\sum_{m=1}^{k}\
\sum_{n_{1}+...+n_{m}=k}\ \sum_{i_{1}\neq...\neq
i_{m}}\binom{k}{n_{1},..,n_{m}}^{*}\hat{x}_{i_{1}}^{n_{1}}...\hat{x}_{i_{m}}^{n_{m}}\sum_{\pi}\hat{y}_{\pi_{i_{1}}}^{n_{1}}...\hat{y}_{\pi_{i_{m}}}^{n_{m}}.$
(17)

The adjustment to the multinomial coefficient is addressed later.

###  3.1 Summation over the Set of Permutations

Beginning with the summation over the permutations, fix $m$,
$i_{1}\neq\cdots\neq 1_{m}$, and $n_{1}+...+n_{m}=k$,

$\sum_{\pi\in\Pi}\hat{y}_{\pi_{i_{1}}}^{n_{1}}...\hat{y}_{\pi_{i_{m}}}^{n_{m}}$
(18)

This sum is independent of the choice of indices $i_{1}\neq...\neq i_{m}$
inherited from the summation over the $x$-coordinates in (17). To see this,
take any two subsets of indices $i_{i}\neq...\neq i_{m}$ and $j_{1}\neq...\neq
j_{m}$ and let $\sigma_{ij}\in\Pi$ be any permutation satisfying
$\sigma_{ij}(i_{l})=j_{l}$ for $l=1,..,m$, ($\sigma_{ij}$ is not unique).
Taking $\Pi$ as a group [13], it follows that $\Pi\circ\sigma_{ij}=\Pi$ and
consequently,

$\displaystyle\sum_{\pi\in\Pi}\hat{y}_{\pi_{i_{1}}}^{n_{1}}...\hat{y}_{\pi_{i_{m}}}^{n_{m}}$
$\displaystyle=$ $\displaystyle=$
$\displaystyle\sum_{\pi\circ\sigma_{ij}\in\Pi}\hat{y}_{(\pi\circ\sigma_{ij})_{i_{1}}}^{n_{1}}\cdots\hat{y}_{(\pi\circ\sigma_{ij})_{i_{m}}}^{n_{m}}$
$\displaystyle=$
$\displaystyle\sum_{\pi\in\Pi}\hat{y}_{\pi_{j_{1}}}^{n_{1}}...\hat{y}_{\pi_{j_{m}}}^{n_{m}}$

Thus we can shift the summation to the first $m$ indices, which we denoted as,

$(n-m)!Y_{(n_{1},...,n_{m})}^{m,k}=\sum_{\pi\in\Pi}\hat{y}_{\pi_{1}}^{n_{1}}...\hat{y}_{\pi_{m}}^{n_{m}}$
(20)

The necessity of $(n-m)!$ will be come clear in a moment. Next, place the
$N=n!$ elements of $\Pi$ into the rows of a matrix and define the operation of
$\Pi$ on the $y$ coordinates as

$\displaystyle\Pi(\hat{y})=\begin{pmatrix}\hat{y}_{\pi_{1}^{1}}&\cdots\cdots&\hat{y}_{\pi_{n}^{1}}\\\
\vdots&\ddots&\vdots\\\
\hat{y}_{\pi_{1}^{N}}&\cdots\cdots&\hat{y}_{\pi_{n}^{N}}\end{pmatrix},$ (21)

(20) is now interpreted as the sum over the rows of the product of the element
in the first $m$ columns raised to the appropriate power. As the order of the
rows does not matter for the summation, rearrange (21) into the canonical
form,

$\Pi^{*}(\hat{y})=\begin{pmatrix}\hat{y_{1}}&\boxed{\begin{minipage}[c][1.85cm][b]{256.0748pt}\centering$$\begin{matrix}\hat{y_{2}}&\boxed{\begin{minipage}[c][1.1cm][b]{227.62204pt}\centering$$\begin{matrix}\hat{y_{3}}&\boxed{\begin{minipage}{202.01474pt}\centering$\cdots$
\@add@centering\end{minipage}}\\\ \vdots\\\ \hat{y_{n}}&\end{matrix}$$
\@add@centering\end{minipage}}\\\ \vdots\\\
\hat{y_{n}}&\boxed{\begin{minipage}[c][0.35cm][b]{227.62204pt}\centering$\vdots$
\ \
\boxed{\begin{minipage}{202.01474pt}\centering$\cdots$\@add@centering\end{minipage}}
\@add@centering\end{minipage}}\end{matrix}$$ \@add@centering\end{minipage}}\\\
\\\
\hat{y_{2}}&\boxed{\begin{minipage}[c][1.6cm][b]{256.0748pt}\centering$$\begin{matrix}\hat{y_{1}}&\boxed{\begin{minipage}[c][0.35cm][b]{227.62204pt}\centering$\vdots$
\ \
\boxed{\begin{minipage}{202.01474pt}\centering$\cdots$\@add@centering\end{minipage}}
\@add@centering\end{minipage}}\\\ \\\
\hat{y_{3}}&\boxed{\begin{minipage}{227.62204pt}\centering$\cdots$
\@add@centering\end{minipage}}\\\ $\vdots$\\\
\hat{y_{n}}&\boxed{\begin{minipage}{227.62204pt}\centering$\cdots$
\@add@centering\end{minipage}}\end{matrix}$$ \@add@centering\end{minipage}}\\\
\vdots\\\ \\\
\hat{y_{n}}&\boxed{\begin{minipage}{256.0748pt}\centering$\cdots$
\@add@centering\end{minipage}}\end{pmatrix}$ (22)

The key observation is that for every distinct sequence of points
$\\{y_{p_{1}},...,y_{p_{j}}\\}$ in the first $j-1$ columns, associated with it
in the $j^{th}$ column are the remaining
$\\{y_{1},...,y_{n}\\}\backslash\\{y_{p_{1}},...,y_{p_{j}}\\}$ each with
multiplicity $(j-1)!$, arising from the remaining $j-1$ columns ignored. Out
of this observation the branching inductive formula is built.

For the case $m=1$, summing along the rows of the canonical form, and a slight
abuse of notation,

$(n-1)!Y_{(k)}^{1,k}=\sum_{\pi\in\Pi}\hat{y}_{\pi_{1}}^{k}=\sum_{Col_{1}}\begin{pmatrix}\hat{y_{1}}^{k}&\boxed{\begin{minipage}{113.81102pt}\centering$\vphantom{\int^{0}}\smash[t]{\vdots}$
\@add@centering\end{minipage}}\\\ \\\
\hat{y_{2}}^{k}&\boxed{\begin{minipage}{113.81102pt}\centering$\vphantom{\int^{0}}\smash[t]{\vdots}$
\@add@centering\end{minipage}}\\\ \vdots\\\
\hat{y_{n}}^{k}&\boxed{\begin{minipage}{113.81102pt}\centering$\vphantom{\int^{0}}\smash[t]{\vdots}$
\@add@centering\end{minipage}}\end{pmatrix}$ (23)

where each factor $\hat{y}_{i}^{k}$ is repeated $(n-1)!$ times, such that

$\displaystyle\sum_{\pi\in\Pi}\hat{y}_{\pi_{1}}^{k}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\hat{y}_{i}^{k}(n-1)!$ (24) $\displaystyle=$
$\displaystyle n\cdot\langle\hat{y}^{k}\rangle(n-1)!$

We now proceed to the branching inductive step. Let $1<m\leq k$ and take non-
zero positive integers $n_{1}+..+n_{m}=k$. We need the following observation,

$\sum_{i\notin\\{i_{1},...,i_{m-1}\\}}\hat{y}_{i}^{n_{m}}\\\
=\Big{[}n\cdot\langle\hat{y}^{n_{m}}\rangle-\hat{y}_{i_{1}}^{n_{m}}-\ldots-\hat{y}_{i_{m-1}}^{n_{m}}\Big{]}.$
(25)

Returning to (20) and reorganizing the sum over distinct sequences of points
in the first $m-1$ columns,

$(n-m)!Y_{(n_{1},...,n_{m})}^{m,k}=\sum_{\pi}\hat{y}_{\pi_{1}}^{n_{1}}\cdots\hat{y}_{\pi_{m}}^{n_{m}},$
(26)

we obtain

$(n-m)!Y_{(n_{1},...,n_{m})}^{m,k}=\sum_{i_{1},..,i_{m-1}}\hat{y}_{i_{1}}^{n_{1}}\cdots\hat{y}_{i_{m-1}}^{n_{m-1}}\sum_{y_{m}\neq
y_{i_{1}},...,y_{i_{m-1}}}\hat{y}_{m}^{n_{m}}(n-m)!$ (27)

which by (25) equals

$=\sum_{i_{1},..,i_{m-1}}\hat{y}_{i_{1}}^{n_{1}}\cdots\hat{y}_{i_{m-1}}^{n_{m-1}}[n\langle\hat{y}^{n_{m}}\rangle-\hat{y}_{i_{1}}^{n_{m}}\cdots-\hat{y}_{i_{m-1}}^{n_{m}}](n-m)!$
(28)

From (28) the branching inductive formula is clear,

$Y_{(n_{1},..,n_{m})}^{m,k}=\ n\langle\hat{y}^{n_{m}}\rangle
Y_{(n_{1},...,n_{m-1})}^{m-1,k-n_{m}}-\sum_{j}Y_{(n_{1}+n_{m}\delta_{ij},...,n_{m-1}+n_{m}\delta_{(m-1)j})}^{m-1,k}$
(29)

The rationale for separating out the factor $(n-m)!$ in (20) is now clear, as
the induction step applies only to the sum over indices and not multiplicity
arsing from the set of permutations. The multiplicity that arises in the
initial sum over permutations is not present in the induction step once the
initial sum is collapsed to a sum over indices alone.

###  3.2 Summation over the $\hat{x}$ terms

Moving to the summation over the $x$-coordinates,

$\sum_{i_{i},..,i_{k}}\hat{x}_{i_{1}}...\hat{x}_{i_{k}}=\sum_{m=1}^{k}\sum_{n_{1}+...+n_{m}=k}\
\ \sum_{i_{1}\neq...\neq
i_{m}}\binom{k}{n_{1},..,n_{m}}^{*}\hat{x}_{i_{1}}^{n_{1}}...\hat{x}_{i_{m}}^{n_{m}},$
(30)

and focusing on

$\sum_{i_{1}\neq...\neq
i_{m}}\hat{x}_{i_{1}}^{n_{1}}...\hat{x}_{i_{m}}^{n_{m}},$ (31)

this is the same exact sum as for the collapsed sum of permutations to a sum
over indices. Hence, starting with the last factor, for $i_{1}\neq\cdots\neq
i_{m-1}$ the sum equals,

$\displaystyle\sum_{i_{1}\neq...\neq
i_{m}}\hat{x}_{i_{1}}^{n_{1}}...\hat{x}_{i_{m}}^{n_{m}}=$ (32)
$\displaystyle=$ $\displaystyle\sum_{i_{1}\neq...\neq
i_{m-1}}\hat{x}_{i_{1}}^{n_{1}}...\hat{x}_{i_{m-1}}^{n_{m-1}}\sum_{i_{m}\neq
i_{1},...,i_{m-1}}\hat{x}_{i_{m}}^{n_{m}}$ $\displaystyle=$
$\displaystyle\sum_{i_{1}\neq...\neq
i_{m-1}}\hat{x}_{i_{1}}^{n_{1}}...\hat{x}_{i_{m}}^{n_{m_{1}}}\Big{[}n\langle\hat{x}^{n_{m}}\rangle-\hat{x}_{i_{1}}^{n_{m}}-...-\hat{x}_{i_{m-1}}^{n_{m}}\Big{]}$

Proceeding in this manner gives us the inductive formula for the $x$ terms,

$X_{(n_{1},..,n_{m})}^{m,k}=\ n\langle\hat{x}^{n_{m}}\rangle
X_{(n_{1},...,n_{m-1})}^{m-1,k-n_{m}}-\sum_{j}X_{(n_{1}+n_{m}\delta_{ij},...,n_{m-1}+n_{m}\delta_{(m-1)j})}^{m-1,k}$
(33)

The adjusted multinomial coefficient represents the number of ways to place
$k$ distinct object into $m$ distinct groups, of respective size
$n_{1},\cdots,n_{m}$, with the caveat that groups of the same size are
indistinguishable. Thus, partition the exponents $(n_{1},\dots,n_{m})$ into
subgroups, $g_{1},\dots,g_{r}$, via the equivalence relation $n_{i}\equiv
n_{j}\iff n_{i}=n_{j}$, and let $d_{i}=|g_{i}|$. The over count that arises,
permuting each group, is equal to $d_{1}!\cdots d_{r}!$, so that

$\binom{k}{n_{1},\dots,n_{m}}^{*}=\
\frac{1}{d_{1}!}\cdots\frac{1}{d_{r}!}\binom{k}{n_{1},\dots,n_{m}}$ (34)

To understand why this adjustment is needed, we proceed with an example. Let
$k=4$, $m=2$, and $n_{1}=n_{2}=2$ and look at the ways the terms in

$\sum_{i\neq j}\hat{x}_{i}^{2}\hat{x}_{j}^{2}$ (35)

arise from the original product of sums on the left side of (30). For fixed
$i\neq j$, We have to pull two $\hat{x}_{i}$ and two $\hat{x}_{j}$ from the
product of the four distinct sums

$\Big{[}\sum_{l}\hat{x}_{l}\Big{]}\Big{[}\sum_{r}\hat{x}_{r}\Big{]}\Big{[}\sum_{s}\hat{x}_{s}\Big{]}\Big{[}\sum_{t}\hat{x}_{t}\Big{]}$
(36)

for which there are $6=\binom{4}{2,2}$ ways to do this

$\displaystyle\hat{x}_{i}\hat{x}_{i}\hat{x}_{j}\hat{x}_{j}\ \ \
\hat{x}_{j}\hat{x}_{j}\hat{x}_{i}\hat{x}_{i}$
$\displaystyle\hat{x}_{i}\hat{x}_{j}\hat{x}_{i}\hat{x}_{j}\ \ \
\hat{x}_{j}\hat{x}_{i}\hat{x}_{j}\hat{x}_{i}$ (37)
$\displaystyle\hat{x}_{i}\hat{x}_{j}\hat{x}_{j}\hat{x}_{i}\ \ \
\hat{x}_{j}\hat{x}_{i}\hat{x}_{i}\hat{x}_{j}$

However, in (37) to form a complete sum each row collapses into a single
complete sum, for a total of $3=\frac{1}{2!}\binom{4}{2,2}$ distinct sums.

Putting everything together gives the desired induction formula for the
moments of Pearson’s distribution over $\Pi(D_{n})$.

∎

For $k=1$, we have

$\langle
r_{\Pi}\rangle=\frac{n^{2}(n-1)!\langle\hat{x}\rangle\langle\hat{y}\rangle}{n\cdot
n!\hat{\sigma}_{x}\hat{\sigma}_{y}},$ (38)

where by definition $\langle\hat{x}\rangle=\langle\hat{y}\rangle=0$ , so that
$\langle r_{\Pi}\rangle=0$ and consequently $\langle r_{\Pi}^{k}\rangle$ are
also central moments.

See Appendix A for the analytic expansions for the first five moments of
$r_{\Pi}$ and numerical validation. Additionally, see Appendix B for a method
organizing the analytic expression for each moment.

##  4 Conclusion

Even without placing these formulas into a statistical framework, it is clear
that a permutation test can not be distribution-free, as some methods for
computing p-value assume. The permutation test depends on the central moments
of the data, as previously investigated in [14]. The exact dependency of the
permutation distribution on data and implications for asymptotic
approximations as sample size increases, can potentially now be investigated
[15, 16, 17]. The problem of approximating a distribution function given a
fixed set of moments has been explored fairly thoroughly[16, 17]. Given a
sufficient number of sample moments for the sample correlation over
permutations of the data, one can estimate the distribution of the sample
correlation over the set of permutations. Such an analysis could potentially
be more efficient with the right implementation, certainly more efficient than
computing permutations of data and their correlation coefficient. We posit
that such a method exists for determining a direct p-value estimate.

## References

  * [1] K. Pearson, “Correlation coefficient,” in Royal Society Proceedings, vol. 58, p. 214, 1895.
  * [2] F. Galton, “Regression towards mediocrity in hereditary stature.,” The Journal of the Anthropological Institute of Great Britain and Ireland, vol. 15, pp. 246–263, 1886.
  * [3] S. M. Stigler, “Francis galton’s account of the invention of correlation,” Statistical Science, pp. 73–79, 1989.
  * [4] J. Lee Rodgers and W. A. Nicewander, “Thirteen ways to look at the correlation coefficient,” The American Statistician, vol. 42, no. 1, pp. 59–66, 1988\.
  * [5] M.-T. Puth, M. Neuhäuser, and G. D. Ruxton, “Effective use of pearson’s product–moment correlation coefficient,” Animal behaviour, vol. 93, pp. 183–189, 2014.
  * [6] E. J. G. Pitman, “Significance tests which may be applied to samples from any populations. ii. the correlation coefficient test,” Supplement to the Journal of the Royal Statistical Society, vol. 4, no. 2, pp. 225–232, 1937.
  * [7] D. S. Chawla, “P-value shake-up proposed,” Nature, vol. 548, pp. 16–17, 2017.
  * [8] M. L. Head, L. Holman, R. Lanfear, A. T. Kahn, and M. D. Jennions, “The extent and consequences of p-hacking in science,” PLoS biology, vol. 13, no. 3, 2015.
  * [9] L. Held and M. Ott, “On p-values and bayes factors,” Annual Review of Statistics and Its Application, 2018.
  * [10] J. I. Krueger and P. R. Heck, “The heuristic value of p in inductive statistical inference,” Frontiers in Psychology, vol. 8, p. 908, 2017.
  * [11] J. I. Krueger and P. R. Heck, “Putting the p-value in its place,” The American Statistician, vol. 73, no. sup1, pp. 122–128, 2019.
  * [12] T. A. Kuffner and S. G. Walker, “Why are p-values controversial?,” The American Statistician, vol. 73, no. 1, pp. 1–3, 2019.
  * [13] M. A. Armstrong, Groups and symmetry.  Springer Science & Business Media, 2013.
  * [14] A. F. Hayes, “Permutation test is not distribution-free,” Psychological Methods, vol. 1, no. 2, p. 184, 1996.
  * [15] N. I. Aheizer and N. Kemmer, The classical moment problem and some related questions in analysis.  Oliver & Boyd Edinburgh, 1965.
  * [16] V. John, I. Angelov, A. Öncül, and D. Thévenin, “Techniques for the reconstruction of a distribution from a finite number of its moments,” Chemical Engineering Science, vol. 62, no. 11, pp. 2890–2904, 2007.
  * [17] R. M. Mnatsakanov, “Hausdorff moment problem: reconstruction of distributions,” Statistics & Probability Letters, vol. 78, no. 12, pp. 1612–1618, 2008.

##  Appendix A Exact Formulas for $k=1,..,5$

From the induction formula every sample moment can be analytically determined.
Below are the exact formulas for the first five sample moments in terms of the
easily computed moments of the data. As in the main body of the paper,
$\hat{\sigma}_{x}$ and $\hat{\sigma}_{y}$ are defined by (13). Additionally,
we employ the following notation for simplicity sake with $\chi_{k}=\langle
x^{k}\rangle$, $\nu_{k}=\langle y^{k}\rangle$, and $\mu_{k,j}=\langle
x^{k}\rangle\langle y^{j}\rangle$, all central moments. The first five sample
moments of the sample correlation are as follows:

$\displaystyle\langle r_{\Pi}^{1}\rangle=$ $\displaystyle\ 0$
$\displaystyle\langle r_{\Pi}^{2}\rangle=$ $\displaystyle\ \frac{1}{(n-1)}$
$\displaystyle\langle r_{\Pi}^{3}\rangle=$ $\displaystyle\
\frac{\mu_{3,3}}{\hat{\sigma}_{x}^{3}\hat{\sigma}_{y}^{3}}\bigg{[}\frac{1}{n^{2}}h_{n,1}+\frac{3}{n^{2}(n-1)}h_{n,2}+\frac{4}{n^{2}(n-1)(n-2)}h_{n,3}\bigg{]}$
$\displaystyle\langle r_{\Pi}^{4}\rangle=$ $\displaystyle\
\frac{1}{\hat{\sigma}_{x}^{4}\hat{\sigma}_{y}^{4}}\bigg{[}\frac{\mu_{4,4}}{n^{3}}h_{n,1}+\frac{4\chi_{4}\nu_{4}}{n^{3}(n-1)}h_{n,2}+\frac{3[n^{2}\hat{\sigma}_{x}^{4}-n\chi_{4}][n^{2}\hat{\sigma}_{y}^{4}-n\nu_{4}]}{n^{5}(n-1)}h_{n,2}$
$\displaystyle+\frac{6[2n\chi_{4}-n^{2}\hat{\sigma}_{x}^{4}][2n\nu_{4}-n^{2}\hat{\sigma}_{y}^{4}]}{n^{5}(n-1)(n-2)}h_{n,3}+\frac{9[2n\chi_{4}-n^{2}\hat{\sigma}_{x}^{4}][2n\nu_{4}-n^{2}\hat{\sigma}_{y}^{4}]}{n^{5}(n-1)(n-2)(n-3)}h_{n,4}\bigg{]}$
$\displaystyle\langle r_{\Pi}^{5}\rangle=$ $\displaystyle\
\frac{1}{\hat{\sigma}_{x}^{5}\hat{\sigma}_{y}^{5}}\bigg{[}\frac{\mu_{5,5}}{n^{4}}h_{n,1}+5\frac{\mu_{5,5}}{n^{4}(n-1)}h_{n,2}+10\frac{[n^{2}\chi_{3}\chi_{2}-n\chi_{5}][n^{2}\nu_{3}\nu_{2}-n\nu_{5}]}{n^{6}(n-1)}h_{n,2}$
$\displaystyle+10\frac{[2n\chi_{5}-n^{2}\chi_{3}\chi_{2}][2n\nu_{5}-n^{2}\nu_{3}\nu_{2}]}{n^{6}(n-1)(n-2)}h_{n,3}$
$\displaystyle+60\frac{[n\chi_{5}-n^{2}\chi_{3}\chi_{2}][n\nu_{5}-n^{2}\nu_{3}\nu_{2}]}{n^{6}(n-1)(n-2)}h_{n,3}$
$\displaystyle+10\frac{[6n\chi_{5}-5n^{2}\chi_{3}\chi_{2}][6n\nu_{5}-5n^{2}\nu_{3}\nu_{2}]}{n^{6}(n-1)(n-2)(n-3)}h_{n,4}$
$\displaystyle+\frac{16[6n\chi_{5}-5n^{2}\chi_{3}\chi_{2}][6n\nu_{5}-5n^{2}\nu_{3}\nu_{2}]}{n^{6}(n-1)(n-2)(n-3)(n-4)}h_{n,5}\bigg{]}$

###  A.1 Validation

In order to demonstrate the validity of these formulas, we randomly generated
datasets of sizes $n=\\{3,\dots,8\\}$ over $100$ trials at each $n$. We
compared the derived moments to the empirical moments computed directly from
$\Pi(D_{n})$, which can be done in the case of small datasets, however for
larger datasets ($n>8$) this becomes computationally impractical.

The mean squared error was computed as:

$MSE=\frac{1}{N_{trials}}\sum\left(\langle r^{k}\rangle_{\Pi(D_{n})}-\langle
r^{k}\rangle_{exact}\right)^{2}$ (39)

This validation procedure was performed using MATLAB with double-precision
floating point computations. The errors tabulated in Table 1 are within
machine epsilon error indicating that the formulas are indeed exact.

| MSE of $k^{th}$ Moment
---|---
Sample Size | 2 | 3 | 4 | 5
$n=3$ | $1.20\mathrm{e}{-32}$ | $2.57\mathrm{e}{-32}$ | $3.99\mathrm{e}{-32}$ | $4.82\mathrm{e}{-32}$
$n=4$ | $6.81\mathrm{e}{-33}$ | $2.46\mathrm{e}{-33}$ | $1.11\mathrm{e}{-32}$ | $3.18\mathrm{e}{-33}$
$n=5$ | $6.06\mathrm{e}{-33}$ | $1.23\mathrm{e}{-33}$ | $5.75\mathrm{e}{-32}$ | $1.55\mathrm{e}{-33}$
$n=6$ | $2.42\mathrm{e}{-32}$ | $8.99\mathrm{e}{-34}$ | $6.00\mathrm{e}{-33}$ | $4.37\mathrm{e}{-34}$
$n=7$ | $1.31\mathrm{e}{-31}$ | $3.50\mathrm{e}{-33}$ | $1.78\mathrm{e}{-32}$ | $1.29\mathrm{e}{-33}$
$n=8$ | $7.38\mathrm{e}{-31}$ | $1.05\mathrm{e}{-32}$ | $5.37\mathrm{e}{-32}$ | $3.17\mathrm{e}{-33}$
Table 1: Validation Error

##  Appendix B Organizing the computation of $\langle r_{\Pi}^{5}\rangle$

Assuming the number of data points $n>5$ and so ignoring the $h_{n,m}$ factor,
from (7) and for $k=5$ we have:

$\displaystyle\langle r_{\Pi}^{5}\rangle$ $\displaystyle=$ $\displaystyle=$
$\displaystyle\frac{1}{n^{k}n!\hat{\sigma}_{x}^{k}\hat{\sigma}_{y}^{k}}\sum_{m=1}^{5}\
\
\sum_{n_{1}+...+n_{m}=5}\binom{5}{n_{1},..,n_{m}}^{*}X_{(n_{1},..,n_{m})}^{m,5}(n-m)!Y_{(n_{1},..,n_{m})}^{m,5}$
(40)

The sum is organized as follows:

For $Z=X$ or $Y$ we have,

$\displaystyle Z_{(1,1,1,1,1)}^{5,5}$ $\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}\hat{z}_{i_{2}}\hat{z}_{i_{3}}\hat{z}_{i_{4}}\hat{z}_{i_{5}}$
$\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}\hat{z}_{i_{2}}\hat{z}_{i_{3}}\hat{z}_{i_{4}}(n\langle
z\rangle-\hat{z}_{i_{1}}-\hat{z}_{i_{2}}-\hat{z}_{i_{3}}-\hat{z}_{i_{4}})$
$\displaystyle=$ $\displaystyle-4Z_{(2,1,1,1)}^{4,5}$ $\displaystyle
Z_{(2,1,1,1)}^{4,1}$ $\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}^{2}\hat{z}_{i_{2}}\hat{z}_{i_{3}}\hat{z}_{i_{4}}$
$\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}^{2}\hat{z}_{i_{2}}\hat{z}_{i_{3}}(n\langle
z\rangle-\hat{z}_{i_{1}}-\hat{z}_{i_{2}}-\hat{z}_{i_{3}})$ $\displaystyle=$
$\displaystyle-Z_{(3,1,1)}^{3,5}-2Z_{(2,2,1)}^{3,5}$ $\displaystyle
Z_{(2,2,1)}^{3,5}$ $\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}^{2}\hat{z}_{i_{2}}^{2}\hat{z}_{i_{3}}$
$\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}^{2}\hat{z}_{i_{2}}^{2}\hat{z}_{i_{3}}(n\langle
z\rangle-\hat{z}_{i_{1}}-\hat{z}_{i_{2}})$ $\displaystyle=$
$\displaystyle-2Z_{(3,2)}^{2,5}$ $\displaystyle Z_{(3,1,1)}^{3,5}$
$\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}^{3}\hat{z}_{i_{2}}\hat{z}_{i_{3}}$
$\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}^{3}\hat{z}_{i_{2}}\hat{z}_{i_{3}}(n\langle
z\rangle-\hat{z}_{i_{1}}-\hat{z}_{i_{2}})$ $\displaystyle=$ $\displaystyle-
Z_{(4,1)}^{2,5}-Z_{(3,2)}^{2,5}$ $\displaystyle Z_{(3,2)}^{2,5}$
$\displaystyle=$ $\displaystyle\sum\hat{z}_{i_{1}}^{3}\hat{z}_{i_{2}}^{2}$
$\displaystyle=$ $\displaystyle\sum\hat{z}_{i_{1}}^{3}(n\langle
z^{2}\rangle-\hat{z}_{i_{1}}^{2})$ $\displaystyle=$ $\displaystyle
n^{2}\langle z^{2}\rangle Z_{(3)}^{1,3}-Z_{(5)}^{1,5}$ $\displaystyle
Z_{(4,1)}^{2,5}$ $\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}^{4}\hat{z}_{i_{2}}$ $\displaystyle=$
$\displaystyle\sum\hat{z}_{i_{1}}^{4}(n\langle z\rangle-\hat{z}_{i_{1}})$
$\displaystyle=$ $\displaystyle-Z_{(5)}^{1,5}$ $\displaystyle Z_{(5)}^{1,5}$
$\displaystyle=$ $\displaystyle\sum\hat{z}_{i_{1}}^{5}$ $\displaystyle=$
$\displaystyle n\langle z^{5}\rangle\text{, where in general
}Z_{(m)}^{1,m}=n\langle z^{m}\rangle$

Combining the terms we have,

$\displaystyle\langle r_{\Pi}^{5}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{n^{k}n!\hat{\sigma}_{x}^{k}\hat{\sigma}_{y}^{k}}\Big{[}(n-1)!X_{(5)}^{1,5}Y_{(5)}^{1,5}$
$\displaystyle+5(n-2)!X_{(4,1)}^{2,5}Y_{(4,1)}^{2,5}+10(n-2)!X_{(3,2)}^{2,5}Y_{(3,2)}^{2,5}$
$\displaystyle+10(n-3)!X_{(3,1,1)}^{3,5}Y_{(3,1,1)}^{3,5}+15(n-3)!X_{(2,2,1)}^{3,5}Y_{(2,2,1)}^{3,5}$
$\displaystyle+10(n-4)!X_{(2,1,1,1)}^{4,5}Y_{(2,1,1,1)}^{4,5}+(n-5)!X_{(1,1,1,1,1)}^{5,5}Y_{(1,1,1,1,1)}^{5,5}\Big{]}$

From which it follows, as seen in Appendix A, that:

$\displaystyle\langle r_{\Pi}^{5}\rangle=$ $\displaystyle\
\frac{1}{\hat{\sigma}_{x}^{5}\hat{\sigma}_{y}^{5}}\bigg{[}\frac{\mu_{5,5}}{n^{4}}+5\frac{\mu_{5,5}}{n^{4}(n-1)}+10\frac{[n^{2}\chi_{3}\chi_{2}-n\chi_{5}][n^{2}\nu_{3}\nu_{2}-n\nu_{5}]}{n^{6}(n-1)}$
$\displaystyle+10\frac{[2n\chi_{5}-n^{2}\chi_{3}\chi_{2}][2n\nu_{5}-n^{2}\nu_{3}\nu_{2}]}{n^{6}(n-1)(n-2)}+60\frac{[n\chi_{5}-n^{2}\chi_{3}\chi_{2}][n\nu_{5}-n^{2}\nu_{3}\nu_{2}]}{n^{6}(n-1)(n-2)}$
$\displaystyle+10\frac{[6n\chi_{5}-5n^{2}\chi_{3}\chi_{2}][6n\nu_{5}-5n^{2}\nu_{3}\nu_{2}]}{n^{6}(n-1)(n-2)(n-3)}$
$\displaystyle+\frac{16[6n\chi_{5}-5n^{2}\chi_{3}\chi_{2}][6n\nu_{5}-5n^{2}\nu_{3}\nu_{2}]}{n^{6}(n-1)(n-2)(n-3)(n-4)}\bigg{]}$


Towards Automated Safety Coverage and Testing for Autonomous Vehicles with
Reinforcement Learning

Hyun Jae Cho, and Madhur Behl

Department of Computer Science

University of Virginia

<EMAIL_ADDRESS>
## Abstract

The kind of closed-loop verification likely to be required for autonomous
vehicle (AV) safety testing is beyond the reach of traditional test
methodologies and discrete verification. Validation puts the autonomous
vehicle system to the test in scenarios or situations that the system would
likely encounter in everyday driving after its release. These scenarios can
either be controlled directly in a physical (closed-course proving ground) or
virtual (simulation of predefined scenarios) environment, or they can arise
spontaneously during operation in the real world (open-road testing or
simulation of randomly generated scenarios).

In AV testing, simulation serves primarily two purposes: to assist the
development of a robust autonomous vehicle and to test and validate the AV
before release. A challenge arises from the sheer number of scenario
variations that can be constructed from each of the above sources due to the
high number of variables involved (most of which are continuous). Even with
continuous variables discretized, the possible number of combinations becomes
practically infeasible to test. To overcome this challenge we propose using
reinforcement learning (RL) to generate failure examples and unexpected
traffic situations for the AV software implementation. Although reinforcement
learning algorithms have achieved notable results in games and some robotic
manipulations, this technique has not been widely scaled up to the more
challenging real world applications like autonomous driving.

The high level contribution of this research project is to build the framework
for automatic scenario testing by using reinforcement learning. In particular,
we use an AV simulator to try to cause a crash for an AV. A safe AV is
expected to drive safely even in rare driving scenarios. We set up two
simulations, one where the RL agent is given a full control over a set of
variables, such as weather, position/velocity of an NPC vehicle, and time of
day, to try to cause a crash for the AV. In the second version, the NPC
vehicle follows a set of fixed waypoints, and the agent is to choose at which
velocity it should be driving at each waypoint to cause a collision.

We discovered that in the first simulation, the reinforcement learning agent
learned that setting a dense fog disabled the AV’s perception module, which
stopped it from moving in the middle of driving. In the second simulation, we
noticed that the agent first let the AV follow its trajectory to its
destination. When it arrived behind the NPC vehicle, the NPC vehicle traveled
at a negative velocity, which the AV did not expect and resulted in a
collision.

During our experiments, we use Baidu Apollo as our AV. We found out that in
the first simulation, the RL agent sets a dense fog that invalidates Apollo’s
perception module from recognizing the traffic lights, making it unable to
drive. In the second experiment, despite Apollo’s robust AD stack, we were
able to cause a direct collision between a non-autonomous NPC vehicle and
Apollo. Therefore, we have revealed potential hazardous driving situations
that current state-of-the-art autonomous vehicles are likely fail to drive
safely. We believe that we have laid the framework for automated safety
testing for autonomous vehicles. We hope that this way of safety testing will
more widely be adopted in the future.

##  1 Background

Today, countless tech corporations and startups are investing billions of
dollars to develop autonomous vehicles (AV). In fact, several businesses will
deploy their vehicles in the next several years [1]. Alphabet’s autonomous
vehicle business, Waymo, recently announced that their commercial self-driving
service, Waymo One, will begin to deploy fully driverless cars [2].

Figure 1: Traditional approaches to maneuvering autonomous vehicles rely on
two broad methods: brute force and formal methods.

There are mainly two approaches in which traditional approaches for training
autonomous vehicles have been taken, depicted in Figure 1. First, autonomous
vehicles are driven as many miles as possible in a brute force way to
encapsulate as many different driving scenarios as possible. However, the
innumerable combinations of the variables involved in driving, such as road
geometry, vehicles, pedestrians, and weather, make it infeasible to encounter
every possible scenario. For example, take Phantom AI’s rear-end collision in
2018 (link to video). Their L2 Advanced Driving Assistance System (ADAS) was
unable to notice a dangerous cut-in by another vehicle and resulted in a
crash. Another popular example is Uber’s fetal crash that killed a pedestrian
in 2018.

When companies advertise the safety of their autonomous vehicles by the number
of miles they have driven, there is one important thing to consider: all
companies drive their vehicles in different environments. An AV that drives
the same high way a million times will not have been exposed to as many unique
driving scenarios as another AV that drove 500 thousand miles in a crowded
urban environment, such as New York City or Downtown Los Angeles. Therefore,
the lack of unified safety testing method for AVs make it ambiguous to compare
safety across multiple autonomous vehicles. For instance, a DMV report in 2018
shows that the number of human engagements in Waymo has decreased dramatically
from 64% in 2015 to 0.09% in 2018, shown in Figure 2. However, not all miles
are equal since every company drives their AVs in different areas. These
examples reveal that the brute force training approach to cover unique
scenarios is just not sufficient.

The other traditional approach for training autonomous vehicles is the use of
formal methods, which reasons about specific driving scenarios and specifies
scenario-specific driving behaviors. However, clearly formal methods is not
capable of generating unique and more complicated scenarios.

Figure 2: Autonomous driving companies such as Waymo drive many miles to
encounter as many unique driving scenarios as possible. In fact, the number of
disengagements has decreased substantially between 2015 and 2018. However, how
much information does number of miles/disengagements tell about safety? Does
this report–which is voluntarily assessed and reported by each
company–necessarily imply that Waymo is safer than its competitors like GM
Cruise or Baidu?

##  2 Motivation

The kind of closed-loop verification likely to be required for AV component
testing is beyond the reach of traditional test methodologies and discrete
verification. Validation puts the autonomous vehicle system to the test in
scenarios or situations that the system would likely encounter in everyday
driving after its release. These scenarios can either be controlled directly
in a physical (closed-course proving ground) or virtual (simulation of pre-
defined scenarios) environment, or they can arise spontaneously during
operation in the real world (open-road testing or simulation of randomly
generated scenarios).

Figure 3: The reinforcement learning agent is tasked with the objective to
maximize the reward associated with causing a collision for a given AV
implementation. In this unique setting, the action space of the agent is the
simulator/scene setup - i.e. the static and dynamic elements of the scene.
After adequate iterations, the agent returns the scenario that is most likely
to cause the AV to fail.

In AV testing, simulation serves primarily two purposes: to assist the
development of a robust autonomous vehicle and to test and validate the AV
before release. A challenge arises from the sheer number of scenario
variations that can be constructed from each of the above sources due to the
high number of variables involved (most of which are continuous). Even with
continuous variables discretized, the possible number of combinations becomes
practically infeasible to test. To overcome this challenge we propose using
reinforcement learning (RL) [3] to generate failure examples and unexpected
traffic situations for the AV software implementation. Although RL algorithms
have achieved notable results in games and some robotic manipulations, this
technique has not been widely scaled up to the more challenging real world
applications like autonomous driving.

The frame of RL is an agent learning through interaction with its environment
driven by an impact (reward) signal. The environment reinforces the agent to
select new actions to improve learning process, hence the name of
reinforcement learning [4]. RL algorithms have achieved notable results in
many domains as games [5, 6] and advanced robotic manipulations [7, 8] beating
human performance. However, standard RL strategies that randomly explore and
learn lose efficiency and become computationally intractable when dealing with
high-dimensional and complex environments [9]. The origin of its difficulty
lies in the important variability inherent to the driving task (e.g.
uncertainty of human behavior, diversity of driving styles, complexity of
scene perception, etc).

In our approach (Figure 3), we rely on the LG Silicon Valley Lab (LGSVL)
simulator [10], a photo-realistic AV simulator that accepts scenario (or
simulation) parameters through a Python API. For an input scenario, the API
provides a way to vary the static (road geometry, position of traffic signs
and lights, number of lanes, etc.) and dynamic parameters (weather, movement
of other agents, pedestrians, etc.) of the simulator. The key insight is that
the parameter space of the simulator, $\Psi$, is the action space for the
reinforcement learning agent. For a given implementation of the AV software
stack, such as Apollo from Baidu [11] or Autoware [12], the agent’s objective
is to choose a set of simulation parameters (static and dynamic) such that it
can cause the AV to fail. For each instance of the simulated scenario, the
agent receives a reward if it can cause the AV to crash. Additionally, it
receives a continuous reward to minimize a safety metric, such as time-to-
collision. The learned policy is a set of simulation parameters that could
cause the AV to fail on the chosen safety metric.

##  3 Related Work

To overcome the shortcomings of the two traditional methodsof AV training,
several recent studies were conducted to train AVs under unexpected, rare
scenarios using reinforcement learning agents and simulation. First, Wang et
al. [13] applied reinforcement learning to train autonomous vehicles for lane
changing tasks in unexpected scenarios. When changing the lane, the lag
vehicle in the target lane can take multiple different actions, such as
yielding by decelerating, accelerating to prevent cut-in, or avoid by changing
lane. To explore the variability, they implemented a model-free reinforcement
learning approach to find the optimal policy when changing lanes. In their
simulation, they includeed some randomness in the parameters, such as the
departure interval and individual speed limits for vehicles. However, their
model was trained on a fixed traffic geometry: a highway segment with fixed
number of lanes, length and width. Therefore, their approach is not scalable
to other road conditions.

In addition, Sarkar et al [14] used the bounded rationality theory to develop
a driving behavior model which generates rate driving scenarios. Their model
used existing driving data to efficiently generate rare driving scenarios that
are strongly correlated to naturalistic driving data. Furthermore, O’Kelley et
al. [15] developed a simulation framework that can test an entire autonomous
driving system. Furthermore, they used adaptive importance-sampling to
evaluate rare-event probability and the probability of an accident in a
simulated driving scenario.

To our knowledge, there has not been a study on using reinforcement learning
and an autonomous vehicle simulator to explicitly target collisions on
existing AVs such as Baidu’s Apollo, which is the goal of this research
project.

##  4 Problem Statement

Let $S$ represent a simulation scenario. In $S$, the RL agent can control a
set of parameters $\Psi$, some of which are dynamic $\Psi_{d}$ and others
static $\Psi_{s}$.

We use RL to generate two simulations: $S_{1}$ and $S_{2}$. In both $S_{1}$
and $S_{2}$, we spawn the AV at the same location and select the same
destination point (Figure 4). The NPC starts from the opposite side of the
crossroad and drives towards the AV. Therefore, the targeted location of
collision is where the trajectories of the two vehicles coincide.

Figure 4: The AV is spawned at the bottom. It follows a trajectory (red) until
it reaches its destination.

The reward function for the RL agent relies on time-to-collision (TTC) between
the NPC vehicle and the Apollo AV.

$TTC(NPC,AV)=\frac{distance(NPC,AV)}{|speed(NPC)-speed(AV)|}$

The denominator is the absolute value of the relative speed of the AV to
NPC’s. As the distance reduces and the two vehicles drive faster toward each
other, TTC will decrease. Therefore, the reward function $R(t)$ for time step
$t$ is defined as

$R(t)=\begin{cases}100&\text{if collision occurs}\\\
\frac{1}{TTC(NPC,AV)}&\text{otherwise}\end{cases}$

This reward function encourages a collision because the agent can gain higher
reward by reducing TTC and by causing a collision. A collision reward of 100
was chosen ineptly.

##  5 Reinforcement Learning Background

Figure 5: Diagram of reinforcement learning [16]. At time step $t$, the agent
takes action $a_{t}$. The environment responds by outputting the new state
$r_{t+1}$ and reward $S_{t+1}$. This cycle continues until the end of training
the model.

Reinforcement learning (RL) is an unsupervised learning algorithm. The
framework in RL involves five main parameters: environment, agent, state,
action, and reward. The environment is the world in which the agent moves. Due
to the unsupervised nature of RL, the agent does not start out knowing the
notion of good or bad actions. At each time step, it performs an action in the
environment that leads to the next state and the corresponding reward (Figure
5). As training progresses, the agent learns which are good and bad actions by
how much reward it can obtain from each action. Therefore, the reward function
plays a vital role for the agent’s ability to learn a policy–a set of
actions–that leads to good performance.

Actor-critic (Figure 6) is a reinforcement learning algorithm that
incorporates two models. The actor model outputs the desired action at each
state. The critic model takes in that action to estimate its value–the
goodness of being in that state. Usually, each model implements a neural
network. At each time step, the weights of the actor model are updated to
perform actions to maximize the value, which the critic model predicts.
Therefore, as training progresses, the actor becomes better at making good
actions and the critic becomes better at estimating the value in each state.
Please refer to pseudo code above for the step by step process.

Figure 6: Architecture of actor critic model [17]. At time step $t$, the agent
takes action $a_{t}$. The environment responds by outputting the new state
$S_{t+1}$ and reward $r_{t+1}$. This cycle continues until the end of training
the model. Actor Critic Algorithm.

⬇

Initialize $V(s)$, the expected reward at state $s$, arbitrarily

and $\pi$ to the policy to be evaluated.

For each episode,

Initialize s.

$a$ $\leftarrow$ action given by actor.

Take action $a$. Observe reward $r$, and next state $s^\prime$.

$V(s) \leftarrow V(s) + \alpha[r + \gammaV(s^\prime)-V(s)]$

$s \leftarrow s^\prime$.

until $s$ is terminal.

##  6 Reinforcement Learning for Safety Coverage Testing

In this section, we describe in detail our process of applying RL to safety
coverage testing. Specifically, we describe the LGSVL simulator, Baidu
Apollo’s AD stack (perception, planning, and control modules), and the
integration of Apollo on the simulator. Finally, we discuss our formulation of
RL to Apollo, which is implemented using the simulator’s Python API.

###  6.1 The LGSVL Simulator

The LGSVL simulator is an open-source, Unity-based autonomous vehicle
simulator. It allows AV developers to test their AD algorithms. Currently, it
supports Autoware and Apollo 3.5 and 5.0 platforms. In addition to the
provided maps and vehicles right out of installation, the simulator allows
creating, editing, and annotating HD maps and vehicles for users.

Figure 7: In our case, we follow the guideline by LGSVL and enable
localization, perception, planning, prediction, control, routing, traffic
light, and transform.

The simulator provides a web user interface (Dreamview) that communicates with
Apollo on simulator via rosbridge. On Dreamview, users can perform a variety
of tasks, including spawning NPC vehicles and pedestrians, visualizing
sensors, and changing weather parameters, including rain, fog, and wetness.
The locations and types of NPC vehicles and pedestrian spawns are chosen
randomly. Using the Python API, NPCs and other parameters can be deterministic
(Figures 8, 9, 10, 11). Users can turn on the modules in Figure 7 to enable
each functionality for the vehicle.

Figure 8: The Python API provides an array of functionalities that can make
desired scenarios. Change and visualize time of day, weather traffic. Figure
9: Visualize sensors and bounding boxes.

Figure 10: Spawn and give waypoints to NPC pedestrians. Figure 11: Spawn NPC
vehicles at desired locations.

We used LGSVL’s BorregasAve map, which renders an actual location in
Sunnyvale, CA. We use the September 2019 release of the simulator.

###  6.2 Baidu Apollo’s AD Stack

Figure 12: Hardware architecture of Baidu’s Apollo. Adequate hardware,
including LiDAR, radar, GPS, and IMU, are equipped to enable the AD stack.
Figure 13: Apollo’s perception pipeline. Apollo uses deep neural networks on
combined data from multiple sensors, such as camera, LiDAR, and radar, for
traffic light and object detection.

Perception, planning, and control are often referred to as the autonomous
driving stack (AD stack). Perception is the module that enables AVs to make
sense of the surroundings. Apollo uses camera images, which are processed and
fed into YOLO v2 [18] object detection algorithm, to locate and classify
pedestrians, traffic signs, lights, other vehicles, and other objects. Then,
it accurately locates itself around the surrounding objects by using a
combination of data from IMU, GPS, radar, and LiDAR, which is then compared to
a previously annotated high definition map (HD map) for increased precision.
Refer to Figure 13 for a detailed pipeline of how Apollo combines data from
multiple sources for detecting traffic lights and objects.

Apollo’s planning module generates a trajectory defined by a series of
waypoints. It calculates the time and velocity it should be driving at when
passing through each waypoint. A curve is then fit on the waypoints to produce
a smooth representation of the trajectory. Planning also observes the
movements of other vehicles and adjusts its waypoints and velocity accordingly
in order to avoid collisions.

The control module is what actually moves the vehicle by controlling
acceleration, braking, and steering angle of Apollo. The waypoints and
corresponding velocities generated from the planning module are given as the
input to this module (Figure 14). Apollo uses both proportional integral
derivative (PID) controller and linear quadratic regulator (LQR) to control
the degree of change in acceleration, braking, and steering angle. PID is
simpler to calculate, but can overshoot the reference trajectory. LQR
complements PID by using the state of the vehicle to minimize lateral error
along the trajectory.

Apollo contains appropriate hardware to support the above modules, shown in
Figure 12.

Figure 14: Diagram showing how the modules from HD maps/localization to CAN
bus work together in Apollo. CAN bus is a communication protocol that allows
components within a vehicle seamlessly. Image source: Udacity’s Self-Driving
Fundamentals: Featuring Apollo course.

###  6.3 Integration of Apollo on the LGSVL Simulator

Apollo receives waypoints from Dreamview via rosbridge and transfers control
commands (accelerate, brake, steering angle) to the simulator which moves the
vehicle. The provided Python API allows sending and receiving websockets to
and from Dreamview to change $\Psi_{d}$, which consists of status of traffic
lights, weather, NPC vehicles, and pedestrians.

###  6.4 RL Formulation to Apollo

We can formulate RL to conduct safety coverage testing of AVs in a simulated
scenario $S$. In every $S$, the Apollo vehicle (EGO) is spawned at a fixed
location and is given a destination with its AD stack on. Therefore, the map
(Borregas Ave), $\Psi_{s}$ (road geometry, position of traffic signs and
lights, number of lanes) and EGO together form the environment of each
scenario. $\Psi_{d}$ form the agent. An agent’s action is a change in
$\Psi_{d}$ at each time step $t$. Finally, the reward function $R(t)$ is set
to the inverse function of time-to-collision (TTC) between the NPC vehicle and
EGO in order to maximize the likelihood of crashing.

The actor of the RL model directs actions for the agent by calculating how
much to change each of $\Psi_{d}$. The critic then estimates the value of each
action by the agent.

##  7 Simulation Results

We have designed our experiment in two versions. First, we allowed the RL
agent to change the following $\Psi_{d}$ in each time step: the position and
velocity of one NPC vehicle, time of day, as well as the amount of rain, fog,
and wetness (continuous values in [0.0, 1.0]). Therefore, the agent’s action
space included changing any combination of these variables. In the second
version, we put restrictions on the agent such that it was only able to change
the velocity of the NPC vehicle. We uniformly sampled a fixed number of
waypoints so that the NPC reached them at each time step. In this version, the
agent’s job was to figure out the velocity to reach each sampled waypoint to
increase the chance to crash with the EGO vehicle. However, the agent was
allowed to travel at negative velocity, in which case it would move backwards.

Figure 15: The RL agent has set a dense fog (fog rate=100%). Figure 16:
Dreamview of the simulator in Figure 15. The EGO vehicle is unable to read the
traffic light, indicated by the ”UNKOWN” sign.

###  7.1 Experiment 1

We discovered that the RL agent will set a dense fog (Figure 15) to invalidate
Apollo’s perception stack from properly perceiving the traffic light (Figure
16), yielding Apollo to stop driving. Although we believe that it is a
logically correct behavior for Apollo to not move, it is also possible that
suddenly stopping in the middle of driving will bring detrimental results of
colliding with other vehicles, which might be moving. This result unveils the
potential inaccuracy of the conventional safety metrics (miles-per-
intervention and the total number of miles a self-driving car has driven)
because we were able to generate a scenario where Apollo’s behavior was not
optimum.

Figure 17: Dreamview of EGO following its trajectory toward its destination.
Note that the green polygon in the path represents the existence of the NPC
vehicle. Figure 18: The NPC vehicle has moved out of the EGO vehicle’s path
and is marked as ”Ignore”.

###  7.2 Experiment 2

When a moving object, such as an NPC vehicle or pedestrian, is on the
trajectory, Apollo will stop and wait until the trajectory is clear of any
obstacle. In Figure 17, Apollo has perceived the existence of an NPC vehicle
(green polygon) on its path. However, as it drove closer, the NPC moved out of
the trajectory. This led Apollo to believe that the NPC was traveling to the
south. Hence, Apollo kept following its trajectory without stopping. However,
the NPC vehicle did not keep driving south; it waited until Apollo was right
behind (Figure 18) and drove at negative velocity, resulting in a crash
(Figures 19, 20).

Figure 19: The NPC vehicle and Apollo collide. Apollo realizes that the NPC
vehicle is driving backwards, indicated by the green arrow, but it was too
late. Figure 20: 19 in the LGSVL simulator view.

##  8 Limitations

The main computational bottleneck in training the RL model came from the LGSVL
simulator. Since Apollo required the camera images to enable its perception
stack, each scenario had to be rendered on the simulator, delaying the time
that could have been spent on exploring more diverse scenarios. Even having
parallel computing would not help because of the rendering problem. The lack
of exclusive command line method of training was the main limitation in this
project.

##  9 Conclusion and Future Work

In this project, we used actor critic reinforcement learning to find edge
cases for Baidu Apollo’s AD stack using the LGSVL simulator. We successfully
generated two edge cases. In the first case, the RL agent controlled a set of
$\Psi_{d}$, such as weather, one NPC vehicle’s positions and velocities, and
time of day to cause a collision between the NPC vehicle and Apollo. Although
it did not specifically crash the two vehicles, the agent impaired the
perception module of Apollo from functioning correctly by creating a dense
fog. Our second edge case, the NPC vehicle was given a set of waypoints, and
agent’s job was to figure out at which velocities it should travel each
waypoint to crash with Apollo. In this scenario, the NPC reached until certain
waypoint, where it stopped and drove backwards when Apollo was behind. The
fact that we were able to generate a collision is significant; it reveals that
the AVs today are potentially vulnerable to unexpected driving scenarios. We
hope that further studies on automated safety testing for autonomous vehicles
can be performed to reveal any other potentially dangerous situations.

Future work can include refining the reward function, adding more NPC vehicles
and pedestrians, and extending the action space of the RL agent.

## References

  * [1] Jon Walker.  The self-driving car timeline – predictions from the top 11 global automakers.  Emerj, Oct 2019.
  * [2] Kirsten Korosec.  Waymo to customers: ’completely driverless waymo cars are on the way’.  TechCrunch, Oct 2019.
  * [3] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [4] Yesmina Jaafra, Jean Luc Laurent, Aline Deruyver, and Mohamed Saber Naceur.  A review of meta-reinforcement learning for deep neural networks architecture search.  arXiv preprint arXiv:1812.07995, 2018.
  * [5] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al.  Human-level control through deep reinforcement learning.  Nature, 518(7540):529, 2015.
  * [6] David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al.  Mastering the game of go with deep neural networks and tree search.  nature, 529(7587):484, 2016.
  * [7] Sergey Levine, Chelsea Finn, Trevor Darrell, and Pieter Abbeel.  End-to-end training of deep visuomotor policies.  The Journal of Machine Learning Research, 17(1):1334–1373, 2016\.
  * [8] Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra.  Continuous control with deep reinforcement learning.  arXiv preprint arXiv:1509.02971, 2015.
  * [9] Niklas Wahlstrom, Thomas B. Schön, and Marc Peter Deisenroth.  From pixels to torques: Policy learning with deep dynamical models.  ArXiv, abs/1502.02251, 2015.
  * [10] Advanced Platform Lab at the LG Electronics America R&D Center (2018), LGSVL Simulator, howpublished = https://github.com/lgsvl/simulator, note = Accessed: 2019-11-20.
  * [11] Haoyang Fan, Fan Zhu, Changchun Liu, Liangliang Zhang, Li Zhuang, Dong Li, Weicheng Zhu, Jiangtao Hu, Hongye Li, and Qi Kong.  Baidu apollo em motion planner.  arXiv preprint arXiv:1807.08048, 2018.
  * [12] Shinpei Kato, Eijiro Takeuchi, Yoshio Ishiguro, Yoshiki Ninomiya, Kazuya Takeda, and Tsuyoshi Hamada.  An open approach to autonomous vehicles.  IEEE Micro, 35(6):60–68, 2015.
  * [13] Pin Wang, Ching-Yao Chan, and Arnaud de La Fortelle.  A reinforcement learning based approach for automated lane change maneuvers.  2018 IEEE Intelligent Vehicles Symposium (IV), Jun 2018.
  * [14] Atrisha Sarkar and Krzysztof Czarnecki.  A behavior driven approach for sampling rare event situations for autonomous vehicles, 2019.
  * [15] Matthew O’Kelly, Aman Sinha, Hongseok Namkoong, John Duchi, and Russ Tedrake.  Scalable end-to-end autonomous vehicle testing via rare-event simulation, 2018.
  * [16] Richard S. Sutton and Andrew G. Barto.  Reinforcement Learning: An Introduction.  The MIT Press, second edition, 2018.
  * [17] Jianing Wei, Anurag Bhardwaj, and wei Di.  Deep Learning Essentials.  Packt, 2018.
  * [18] Joseph Redmon and Ali Farhadi.  YOLO9000: better, faster, stronger.  CoRR, abs/1612.08242, 2016.