# Modelling Periodic Measurement Data
Having a Piecewise Polynomial Trend Using the
Method of Variable Projection

Johannes Handler, Dimitar Ninevski and Paul O’Leary This work was partially
funded by the COMET program within the K2 Center “Integrated Computational
Material, Process and Product Engineering (IC-MPPE)” (Project No 859480). This
program is supported by the Austrian Federal Ministries for Transport,
Innovation and Technology (BMVIT) and for Digital and Economic Affairs (BMDW),
represented by the Austrian research funding association (FFG), and the
federal states of Styria, Upper Austria and Tyrol.Johannes Handler, Dimitar
Ninevski and Paul O’Leary are with the Chair of Automation, Department of
Product Engineering, University of Leoben, A8700 Leoben, Austria (email:
<EMAIL_ADDRESS>
###### Abstract

This paper presents a new method for modelling periodic signals having an
aperiodic trend, using the method of variable projection. It is a major
extension to the IEEE-standard 1057 by permitting the background to be time
varying; additionally, any number of harmonics of the periodic portion can be
modelled. This paper focuses on using B-Splines to implement a piecewise
polynomial model for the aperiodic portion of the signal. A thorough algebraic
derivation of the method is presented, as well as a comparison to using global
polynomial approximation. It is proven that B-Splines work better for
modelling a more complicated aperiodic portion when compared to higher order
polynomials. Furthermore, the piecewise polynomial model is capable of
modelling the local signal variations produced by the interaction of a control
system with a process in industrial applications. The method of variable
projection reduces the problem to a one-dimensional nonlinear optimization,
combined with a linear least-squares computation. An added benefit of using
the method of variable projection is the possibility to calculate the
covariances of the linear coefficients of the model, enabling the calculation
of confidence and prediction intervals. The method is tested on both real
measurement data acquired in industrial processes, as well as synthetic data.
The method shows promising results for the precise characterization of
periodic signals embedded in highly complex aperiodic backgrounds. Finally,
snippets of the m-code are provided, together with a toolbox for B-Splines,
which permit the implementation of the complete computation.

###### Index Terms:

B-Splines, data modelling, IEEE 1057, signal separation, variable projection
method.

##  I Introduction

A portion of this paper was originally published at the I2MTC conference in
2021 [1]. There, the authors showed the advantages of applying the _method of
variable projection_ [2, 3] to identifying the four-parameter model for a sine
wave, see the IEEE-standard 1057 [4] for an exact definition of the model.
This model is commonly used in waveform digitizers, e.g. digital
oscilloscopes, analog to digital converters etc [5, 6]. However, the authors
previous work was primarily motivated by the widespread availability of low-
cost MEMS accelerometers and gyroscopes [7]; which, commonly create such
signals. The ubiquitous applicability of these sensors [8, 9, 10] calls forth
the need for stable and computationally efficient methods.

The four-parameter model assumes a constant background term $d$. This paper
extends the method to a much wider class of background signals; that is, the
background $d(t)$ is a function of time,

$y_{m}(t)=a\,\sin(\omega t+\phi)+d(t).$ (1)

The aim is to make the new methods much more widely applicable in the
instrumentation of industrial processes.

A periodic signal, without an integer number of cycles in the measurement
period $T$, combined with an aperiodic background $d(t)$, is also the
fundamental issue behind spectral leakage when computing Fourier spectra [11].
Consequently, the methods presented here are, not only relevant to industrial
instrumentation, but also relevant to many cases of signal processing.
Traditionally, this issue has been addressed using windowing techniques, see
the seminal paper by Harris [12].

###  I-A Difficulty of separation

Figure 1: Two real sensor signals measured in different industrial processes.
Note there is a global periodic portion to the signal on top of an aperiodic
background. The aperiodic portion is due to the process control making changes
to the actuators driving the process.

In Figure 1, two real sensor signals measured in different industrial
processes111For confidentiality reasons we are not permitted to name the
process behind these measurements. are shown. Both exhibit a combination of
periodic and aperiodic portions. The measurement data is typical for many
industrial applications where there is a control process in the background
making changes to the actuators driving the process. The goal here is to
separate these two components reliably and characterize the periodic portion
exactly.

The magnitude of the FFT for the two sensor can be seen in Figure 2. Note: The
frequency range has been cropped to show the region where the signals are to
be observed. The FFT is highly attractive due to its
$\mathcal{C}(n)=n\,log(n)$ computational efficiency [13]. However, as can be
seen in Figure 2, the spectral leakage from the non-integer number of cycles
and from the background signal make the reliable characterization of the
periodic portion intractable. Additionally, the sharp changes in the signal,
due to the control process, may lead to a Gibbs error with periodic portions.

Figure 2: The magnitude of the FFT for the two signals shown in Figure 1. The
frequency range has been cropped to show the region where the signals are to
be observed and the DC component removed. The respective sampling frequencies
$f_{s}$ and measurement periods $T$ are also shown.

The problem remains the identification of a suitable model for the $d(t)$ term
in Equation 1 and then the implementation of a stable, efficient and reliable
method to determine the coefficients that lead to a good approximation of the
signal.

Similar problems can be found in literature: that is, problems where the model
equation is a linear combination of nonlinear functions appear frequently and
are commonly solved using the variable projection method. Such problems can be
found in power system modal identification [14], GPS positioning [15, 16],
computer vision [17, 18, 19], calibration of measurement data for
electromagnetic data inversion [20], etc. Many of these problems have a
combination of exponential and trigonometric functions as nonlinear basis
functions, such as [14, 21]. Others [16, 15, 22] have looked at the case when
there are some linear functions, e.g. polynomials, in the nonlinear basis
matrix and how to effectively use this to further simplify the calculation. In
[22], the authors use a combination of polynomials and trigonometric functions
to separate the periodic and trend portions of the signal; however, issues
arise when the background signal $d(t)$ becomes more complicated and higher
degree polynomials are required. Later in this paper it is proven that higher
order polynomials, for the background, also start to model portions of the
periodic function: This can be seen from the Taylor expansion. Consequently,
higher order polynomials can lead to a less accurate characterization of the
frequency and magnitude of the periodic portion.

The main contributions of this paper are:

  1. 1.

To propose a piecewise polynomial as a model for the background signal $d(t)$.
This model is compatible and consistent with industrial measurements where the
process is being modified by the intervention of actuators. The process
control is locally, but not globally, piecewise stationary. This property can
be utilized to obtain a better severation from the global periodic oscillation
we wish to characterize.

  2. 2.

The B-Spline implementation of the piecewise polynomials, permits a direct
combination with periodic basis functions; in this manner the method of
variable projection can be applied. This splits the originally $n$-dimensional
optimization into a one-dimensional nonlinear minimization to determine the
frequency $\omega$, followed by an $n-1$ dimensional linear problem given
$\omega$. For the linear portion the methods, by definition, ensure that the
global minimum given $\omega$ is found. Additionally, the covariance of the
linear coefficients is simply computed.

  3. 3.

A thorough analysis of the theoretical background to the models used and the
methods applied to obtain the model coefficients is presented.

  4. 4.

Snippets of m-code are provided which implement the proposed methods. An
m-code library for the implementation of the B-Spline bases is provided.

  5. 5.

Results are presented which validate the method using synthetic data sets with
known properties.

  6. 6.

The results from the successful application of the new approach to the
analysis of industrial measurement data are presented, from two differing
processes.

##  II Signal models

The goal now is to extend the four-parameter model [4] to include a time
dependent background $d(t)$. This will accommodate more complex signals, while
maintaining the ability to accurately characterize the periodic portion.

###  II-A Piecewise polynomial background

In this paper we propose a piecewise polynomial [23] for the background signal
$d(t)$, implemented via a B-Spline model [24, Chapter 2]. The justification
for this is: the control signals, which are causal for the aperiodic portion,
are in general piecewise $C^{0}$ continuous. However, due to the dynamics of
the process, the observed sensor signals tend to be $C^{1}$ or $C^{2}$
continuous. It is the differential equation describing the dynamics of the
process, that govern the order $n$ of the continuity $C^{n}$ observed in the
sensor signals. The piecewise $C^{n}$ continuity of the background signal
ensure that piecewise polynomials are suitable models for $d(t)$. A further
advantage in industrial measurements is that the break points, required for
the piecewise polynomials, are known from the control signals.

The B-Spline implementation has been chosen, over the pp-spline, since it
permits the direct combination of the basis functions $\bm{B}_{s}$ for the
spline and $\bm{B}_{p}$ for the periodic portions, into one consistent set of
matrix equations. In this manner the method of variable projection [2, 3] can
be applied, during the optimization, to separate the nonlinear from the linear
portion of the computation. This improves the stability of the optimization
process and reduces the computational complexity. The possibility of a
piecewise continuous background $d(t)$, extends greatly the areas where this
methods can be applied in the analysis of measurement signals. The four
parameter model [4], considered in the previous paper [1], is now a special
case also covered by the model proposed here.

###  II-B B-Spline notation

To compute a complete set of spline basis functions we need to know the vector
of points $\bm{x}$, where the spline is to be evaluated, the desired degree of
the spline $\delta$ and a vector of knots $\bm{\kappa}$. Consequently, the
notation222This notation is compatible with the B-Spline m-code implementation
made available at MATLAB-FileExchange [25]. for the spline should be of the
form $\bm{B_{s}}(\bm{x},\delta,\bm{\kappa})$. The spline portion of the model
$\bm{y}_{s}$ is now computed as a linear combination of these basis functions
with the coefficient vector $\bm{\beta}$, i.e.,

$\bm{y}_{s}=\bm{B_{s}}(\bm{x},\delta,\bm{\kappa})\,\,\bm{\beta}.$ (2)

The number of basis functions required and with this the length of the
coefficient vector $\dim(\beta)$ is dependent of $\delta$ and the number of
knots, i.e.,

$\dim(\beta)=\dim(\bm{\kappa})+\delta-1$ (3)

###  II-C Periodic notation

The periodic portion is modelled as a base frequency $\omega$, together with a
number of harmonics $\nu$ that need to be considered. The use of harmonics
permits the modelling of more complex periodic signals. The complete notation
required for the periodic basis functions is $\bm{B_{p}}(\bm{x},\omega,\nu)$.
The coefficient vector $\alpha$ has two components per harmonic,

$\dim(\bm{\alpha})=2\,\nu,$ (4)

since one sine and one cosine is required per harmonic. The periodic portion
of the signal is computed as,

$\bm{y}_{p}=\bm{B_{p}}(\bm{x},\omega,\nu)\,\,\bm{\alpha}.$ (5)

###  II-D Complete signal model

The complete signal model $\bm{y}_{m}=\bm{y}_{p}+\bm{y}_{s}$, which in matrix
form can now be written as,

$\bm{y}_{m}=\begin{bmatrix}\bm{B_{p}}(\bm{x},\omega,\nu)&\bm{B_{s}}(\bm{x},\delta,\bm{\kappa})\end{bmatrix}\,\,\begin{bmatrix}\bm{\alpha}\\\
\bm{\beta}\end{bmatrix}.$ (6)

To obtain a less cumbersome notation we shall define some simplified
notations: since the locations of the samples $\bm{x}$, $\nu$ the number of
harmonics, $\delta$ the degree of the polynomial spline and $\bm{\kappa}$ the
vector of break points do not change during the optimization, we shall define

$\bm{B}(\omega)\triangleq\begin{bmatrix}\bm{B_{p}}(\bm{x},\omega,\nu)&\bm{B_{s}}(\bm{x},\delta,\bm{\kappa})\end{bmatrix}\,\,\text{given}\,\,\bm{x},\nu,\delta,\bm{\kappa},$
(7)

indicating that $\bm{B}(\omega)$ is a matrix of basis functions nonlinear in
$\omega$ and $\omega$ is relaxant during the optimization. Additionally,

$\bm{\gamma}\triangleq\begin{bmatrix}\bm{\alpha}\\\ \bm{\beta}\end{bmatrix}.$
(8)

Now, the notation for the model calculation simplifies to,

$\bm{y}_{m}=\bm{B}(\omega)\,\,\bm{\gamma}.$ (9)

The least squares approximation of the signal $\bm{y}$ by the model
$\bm{y}_{m}$, now requires us to minimize,

$\min_{\omega,\bm{\gamma}}\|\bm{y}-\bm{B}(\omega)\,\bm{\gamma}\|^{2}_{2}.$
(10)

Attacking this task directly leads to a very high dimensional nonlinear
optimization problem, with all the associated numerical difficulties.
Fortunately, the problem separates into a linear combination of basis
functions which are nonlinear in $\omega$. This makes the system of equations
suitable for the application of variable projections.

##  III Method of Variable Projection

The method of variable projection [2, 3] is a method for solving separable
nonlinear least squares problems. It is characterized by the model being a
linear combination of nonlinear functions, meaning

$y_{m}=\alpha_{1}b_{1}\left(\bm{\delta}\right)+\alpha_{2}b_{2}\left(\bm{\delta}\right)+\ldots+\alpha_{k}b_{k}\left(\bm{\delta}\right),$
(11)

where $\bm{\delta}$ is a vector of $j$ nonlinear coefficients,
$b_{i}\left(\bm{\delta}\right)$ are the basis functions, which are nonlinear
in $\bm{\delta}$ and $\alpha_{i}$ are the $k$ linear coefficients.
Concatenating the individual basis functions into the matrix of basis
functions, yields

$\bm{B}(\bm{\delta})=[\bm{b}_{1}(\bm{\delta}),\bm{b}_{2}(\bm{\delta}),\ldots,\bm{b}_{k}(\bm{\delta})].$
(12)

The notation $\bm{B}(\bm{\delta})$ indicates that the contents of the matrix
$\bm{B}$ are dependent on $\bm{\delta}$. Now defining the coefficient vector

$\bm{\alpha}\triangleq[\alpha_{1},\alpha_{2},\ldots\alpha_{k}]^{\mathrm{T}}$
(13)

leads directly to the matrix vector equation,

$\bm{y}_{m}=\bm{B}(\bm{\delta})\,\bm{\alpha}.$ (14)

Given an estimate for $\bm{\delta}$, which yields a stationary value for the
cost function, a least squares estimate for $\bm{\alpha}$ is obtained from,

$\bm{\alpha}=\bm{B}^{+}(\bm{\delta})\,\bm{y},$ (15)

whereby, $\bm{B}^{+}(\bm{\delta})$ denotes the Moore-Penrose pseudo inverse of
$\bm{B}(\bm{\delta})$. Now substituting Equation 15 into 14 one obtains,

$\bm{y}_{m}=\bm{B}(\bm{\delta})\,\bm{B}^{+}(\bm{\delta})\,\bm{y}.$ (16)

Defining
$\bm{P}(\bm{\delta})\triangleq\bm{B}(\bm{\delta})\,\bm{B}^{+}(\bm{\delta})$,
this is the projection onto the subspace spanned by the basis functions
contained in $\bm{B}(\bm{\delta})$. This projection varies with $\bm{\delta}$,
hence the name, method of variable projection. The residual vector is computed
as

$\bm{r}=\bm{y}-\bm{y}_{m}$ (17)

leading to the cost function $E(\bm{\delta})$, defined as the sum of squares
of the residual $\bm{r}$, to be calculated as

$\displaystyle E(\bm{\delta})$ $\displaystyle=\|\bm{r}\|^{2}$ (18)
$\displaystyle=\|\bm{y}-\bm{B}(\bm{\delta})\,\bm{B}^{+}(\bm{\delta})\,\bm{y}\|^{2},$
(19)
$\displaystyle=\|\left\\{\bm{I}-\bm{B}(\bm{\delta})\,\bm{B}^{+}(\bm{\delta})\right\\}\,\bm{y}\|^{2}.$
(20)

This is called the _variable projection functional_ (VPF). Note that the
estimation of $\bm{\delta}$ is now a nonlinear least squares problem in $j$
dimensions; whereas, the original problem has $j+k$ coefficients. Equation 18
permits an explicit calculation of the cost function $E(\bm{\delta})$ as a
function of $\bm{\delta}$ for a given measurement $\bm{y}$. In Section 4 of
[2], Golub provides the formal proofs required to determine that the Fréchet
derivatives over the pseudo-inverse and projection yield the gradient of the
cost function. This implies that gradient based nonlinear solvers can be used
to find the value of $\bm{\delta}$ which minimizes the cost function
$E(\bm{\delta})$.

This method is particularly advantageous for the case at hand, since the cost
function is only nonlinear in the variable $\omega$. From Equation 10 we
obtain the variable projection (cost) function,

$E(\omega)=\|\bm{y}-\bm{B}(\omega)\,\bm{B}^{+}(\omega)\,\bm{y}\|^{2}_{2}.$
(21)

That is, the problem has reduced to a one-dimensional nonlinear optimization
problem.

##  IV Covariance propagation

A further advantage of the method of variable projection is that it yields a
direct method of calculating the covariance propagation for the linear
coefficients. Given $\bm{\gamma}=\bm{B}^{+}(\omega)\,\bm{y}$, a linear
mapping, then $\bm{\Lambda}_{\gamma}$, the covariance matrix of $\bm{\gamma}$,
can be computed as [26],

$\bm{\Lambda}_{\gamma}=\bm{B}(\omega)^{+}\,\bm{\Lambda}_{y}\,{\bm{B}(\omega)^{+}}^{\mathrm{T}}.$
(22)

whereby, $\bm{\Lambda}_{y}$ is the covariance matrix of the data vector
$\bm{y}$. If $\bm{y}$ is perturbed by i.i.d. Gaussian noise with the standard
deviation $\sigma_{y}$, then the equation becomes

$\bm{\Lambda}_{\gamma}=\sigma^{2}_{y}\,\bm{B}(\omega)^{+}\,{\bm{B}(\omega)^{+}}^{\mathrm{T}}.$
(23)

Assuming the model is bias free333In applications the Gaussian nature of the
residual vector $\bm{r}$ should be verified to ensure that this assumption is
valid., an estimate for $\sigma_{y}^{2}$ can be computed from the residual
vector $\bm{r}$ as follows:

$\sigma_{y}^{2}=\frac{1}{n-n_{df}}\,\|\bm{r}\|_{2}^{2},$ (24)

where by $n_{df}$ denotes the number of degrees of freedom. In this manner,
the covariance of the linear coefficients can be computed directly.

##  V Interaction between the portions

The question now is: what would be an ideal signal for $d(t)$? Ideally the
model for the background $d(t)$ would be fully orthogonal to the periodic
portion, to ensure there is no interaction between the components. The inner
product of two vectors
$c=\bm{a}^{\textrm{T}}\,\bm{b}=|\bm{a}|\,|\bm{b}|\,\cos(\phi)$, where $\phi$
is the angle between the two vectors, i.e., it is a measure for their
orthogonality. The inner product for a continuous $d(t)$ wrt. to a single
periodic component $y_{p}(t)=\mathrm{e}^{-j\omega t}$ can be calculated as,

$c=\int_{-\infty}^{\infty}d(t)\,\mathrm{e}^{-j\omega t}\,\mathrm{d}t.$ (25)

This corresponds to the fourier coefficients for $d(t)$ up to a factor of
scale. Consequently, for discrete $d(t)$ we can apply the FFT to obtain a
numerically efficient estimate for the interaction between the two portions of
the signal.

The Weierstrass theorem [27] states that, in a finite interval, any function
can be approximated by a polynomial with the desired accuracy, given a
polynomial of sufficient degree. This is also true for periodic signals;
consider the Maclaurin series for $\sin(x)$ and $\cos(x)$,

$\displaystyle\sin\left(x\right)$
$\displaystyle=x-\frac{x^{3}}{3!}+\frac{x^{5}}{5!}-...=\sum_{k=0}^{\infty}\left(-1\right)^{k}\frac{x^{2k+1}}{\left(2k+1\right)!}$
(26) $\displaystyle\cos\left(x\right)$
$\displaystyle=1-\frac{x^{2}}{2!}+\frac{x^{4}}{4!}-...=\sum_{k=0}^{\infty}\left(-1\right)^{k}\frac{x^{2k}}{\left(2k\right)!}.$
(27)

Consequently, a global polynomial, as a model for $d(t)$, will have an
interaction with the periodic portion. The level of interaction can be
computed by creating the Vandermonde matrix $\bm{V}$ for $-1\leq x\leq 1$, and
applying the FFT to its columns. This yields the spectra of the monomials, see
Figure 3.

Figure 3: Top: Fourier spectra of the monomials $x^{0}\ldots x^{5}$ in the
interval $-1\leq x\leq 1$. Bottom: The Fourier spectra for the $n=6$ B-Spline
basis functions used for the synthetic data, see Section VI-A. Note: there is
significantly less interaction of the B-Splines with the periodic signal.

In the case of B-Splines[28] the interaction depends on the locations of the
knots. Note that when using splines, determining the locations of the
breakpoints (knots) is very important and often very challenging. Here,
however, the focus is on the variable projection method and for further
information on the placement of knots, the reader is referred to [29, 30].
Fortunately, in the application fields driving this work, the locations of the
knots are known from the control system data.

At this point the knots associated with the synthetic data, see Section VI-A,
are used; since this permits a comparison of B-Spline bases with monomials.
The interaction between the B-Splines and the periodic signal is also shown in
Figure 3.

The Fourier spectra of the monomials and the B-Spline bases, show that there
is significantly less interaction between the periodic signal and the
B-Splines. Furthermore, the interaction will diminish with increasing number
of knots, since the B-Spline become more local in nature. This is the reason
why using polynomials, of higher degree, to model a trend globally does not
always provide good results [22]. Consequently, the B-Splines can be
considered as the preferred solution for $d(t)$ over global polynomials.

##  VI Numerical testing

The proposed method was tested on different datasets: a synthetic dataset, and
two real sensor signals emanating from an industrial process, where a control
process is making changes to the actuators driving the process.

###  VI-A Synthetic examples

For the purposes of testing, a synthetic data set was generated as the sum of
a spline and a trigonometric function, as shown in Figure 4. Here the knot
locations $\kappa$ were known a-priori. The piecewise polynomial trend is
defined by five knots. The degree of the spline is chosen to be $\delta=2$.
The periodic component consists a sine with frequency $\omega=36.96$ and its
first two harmonics. Finally, the synthetic signal is the sum of those two
components superimposed with i.i.d. Gaussian noise with a standard deviation
of $\sigma=0.05$ (Figure 4, bottom).

Figure 4: The two components of the synthetic test signal and the final test
signal defined for a normalized time vector $x$: (top) the piecewise
polynomial trend, defined by 5 knots and with the degree of the polynomial
spline being $\delta=2$, (middle) the periodic component consisting of a sine
wave with frequency $\omega=36.96$ and its first two harmonics, (bottom) the
sum of those two components superimposed with i.i.d. Gaussian noise with a
standard deviation of $\sigma=0.05$.

In general, the result of a numerical optimization strongly depends on the
initial value. Here, the FFT provides a good starting point for the numerical
iteration procedure. However, applying the FFT directly on the measurement
signal $y$ does not provide good results; since the aperiodic background
$d(t)$ leads to a major distortion of the spectrum, see Figure 6 (top). In
order to reduce this influence the signal $y$ is first approximated solely by
the B-Spline basis functions,

$\bm{y}_{s}=\bm{B_{s}}\,\bm{B_{s}}^{+}\,\bm{y}.$ (28)

The residual of this approximation, $\bm{r}_{s}=\bm{y}-\bm{y}_{s}$, contains
primary periodic portions, as shown in Figure 5. The frequency corresponding
to the maximum magnitude in the spectrum of the residual $r_{s}$ is the base
frequency we are looking for, subject to Gibbs error [11]. However, this value
is a good initial value for the nonlinear optimization process to get a more
accurate approximation for the base frequency.
In the presented example this gives an initial frequency of
$\omega_{init}=37.699$ and lead to an estimation of $\hat{\omega}=36.884$.

Figure 5: For the computation of a suitable initial value, the measurement
data is solely approximated by the B-Spline basis functions: (top) synthetic
data and the approximating B-Spline, (bottom) the residual of the
approximation. Figure 6: The magnitude of the FFT for the two signals shown in
Figure 5: (top) The spectrum of the synthetic data, (bottom) the spectrum of
the residual.

The final result, as well as the residual, is shown in Figure 7. These results
demonstrate the ability of the computational method to successfully separate
the components of the signal.

Figure 7: The result of the data approximation by the method of variable
projection: (top) The noisy data, and the estimated function, (middle two) the
resulting spline and periodic components, (bottom) the residual.

Additionally, for this example, the covariances of the linear coefficients
were calculated. Since a spline of degree $2$ with $5$ knots was used,
according to Equation 3 this results in $6$ spline coefficients denoted by
$\beta_{1},\beta_{2},\ldots,\beta_{6}$. The covariances of those coefficients
can be seen in Table I. Note that this model has more linear coefficients,
namely the linear coefficients in front of the trigonometric functions,
however these were not displayed here for conciseness.

Cov | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\beta_{5}$ | $\beta_{6}$
---|---|---|---|---|---|---
$\beta_{1}$ | 0.132 | -0.053 | 0.023 | -0.010 | 0.009 | -0.007
$\beta_{2}$ | -0.053 | 0.091 | -0.044 | 0.019 | -0.015 | 0.008
$\beta_{3}$ | 0.023 | -0.044 | 0.063 | -0.029 | 0.022 | -0.010
$\beta_{4}$ | -0.010 | 0.019 | -0.030 | 0.053 | -0.041 | 0.019
$\beta_{5}$ | 0.009 | -0.015 | 0.022 | -0.041 | 0.100 | -0.054
$\beta_{6}$ | -0.007 | 0.008 | -0.010 | 0.019 | -0.054 | 0.105
TABLE I: Table of covariances for the test case shown in Figure 7. The results
are scaled by $10^{3}$ since they were very small.

###  VI-B Industrial measurement data

The presented algorithm was further tested on data sets emanating from
different industrial processes.444Due to confidentiality reasons the original
data sets were anonymized. A coordinate transformation is applied to $x$ to
yield dimensionless values. The x-scaling has been chosen using apriory
knowledge, so that the expected periodicity has a period of $T_{e}=1$. As
mentioned before, choosing the right locations of the breakpoints $\kappa$ and
the degree of the spline $\delta$ has a strong influence on the quality of the
result. Here, the knots are chosen to correspond to time points when the
control system influences the process via an actuator. This has the advantage
that the knot locations can be derived directly from the control signals.

In the first example, shown in Figure 8, it can be seen that the aperiodic
component $y_{s}$ exhibits an oscillatory behavior at around $x=8$. This is
due to the local activation of the process dynamics by the process control.
Due to the local nature of splines they are well suited to model this
behavior. The natural frequency of the process dynamics is close to the
frequency of the global periodic component $y_{p}$. Despite this fact the
proposed algorithm is still in the position to successfully separate the two
signal components. Also the period being identified as $T_{i}=1.04$ is in the
expected range. This clearly speaks for the robustness of the method.

Data from a second process is presented in Figure 9. Here the resulting period
after optimization is $T_{i}=0.996$. In this measurement data are some areas,
e.g. $x=28$ and $x=56$ where the residual is relatively large; this requires
further discussion. Considering the residual in more detail: the disturbances
are primarily impulse in nature and not periodic; this indicates that $d(t)$
may need to be modified to improve the modelling. The local nature of
B-Splines permits the insertion of additional knots [28]. The new breakpoints
should be inserted where the local integral of $r(t)^{2}$ is a maximum: This
is an area where further research is justified. Furthermore, this would open
the possibility to identify and quantify extraneous influences not emanating
from the control.

Figure 8: Measurement data from an industrial process. (top) Measurement
signal and the approximated function, whereby the period was identified to be
$1.04$. (middle two) The spline and periodic portion of the signal obtained
using the method of variable projection. (bottom) Residual of the
approximation.  Figure 9: Another industrial measurement data set. (top)
Measurement signal and the approximated function, whereby the cycle duration
was identified to be $0.996$. (middle two) The spline and periodic portion of
the signal obtained using the method of variable projection. (bottom) Residual
of the approximation.

##  VII Code implementation

Here the salient snippets of code required to implement the method of variable
projection for modelling a periodic signal having a piecewise polynomial
background are presented. The necessary support functions555Harmonic basis
functions are choosen instead of a single sine wave to model more complex
periodic signals. harmonicBasisand bSplineBasis are made available at MATLAB-
FileExchange [25]. The first piece of m-code is a function that is used to
compute the cost function $E(\omega)$ for the method of variable projection,
see Listing 1.

⬇

1function cost = varproCost( omega, x, y, nu, delta, kappa )

2%

3% Setup the basis functions

4Bp = harmonicBasis( x, omega, nu );

5Bs = bSplineBasis( x, delta, kappa );

6B = [ Bp, Bs ];

7% Compute the projection

8ym = B * ( B \ y );

9% Evaluate the residual

10r = y - ym;

11% Compute the corresponding cost

12cost = norm( r );

Listing 1: Variable projection cost function

As in [31], here a standard high quality generic nonlinear iteration
process666This is a common approach, since such solvers are available in most
numerical computation packages. is availed to perform the optimization. The
code required for this portion is shown in Listing 2. An autonomous function
fun is defined with one calling parameter, omega and one return parameter
cost. Additionally, there are five passively passed parameters x, y, nu, delta
and kappa: the time vector and measurement vector, respectively, as well as
the parameters defining the B-Spline bases and the number of harmonics.

⬇

1% Define the anonymous for the cost

2fun = @(omega) varproCost( omega, x, y, nu, delta, kappa );

3% Wrap with the nonlinear iteration.

4omega = lsqnonlin( fun, omega_Init );

Listing 2: Wrapping with a nonlinear solver

At the end of running the code in Listing 2, the optimal value for $\omega$ is
obtained. The remaining linear coefficients $\bm{\gamma}$, since $\omega$ is
available, can be computed according to Equation 15, this is implemented in
Listing 3.

⬇

1Bp = harmonicBasis( x, omega, nu );

2Bs = bSplineBasis( x, delta, kappa );

3B = [ Bp, Bs ];

4% Compute the projection

5cfs = B \ y ;

Listing 3: Compute the linear parameters.

Finally, if required the covariances of the linear coefficients can be
computed according to the m-code in Listing 4.

⬇

1 nrParams = size( B, 2 ) + 1;

2 % Number of degrees of freedom

3 df = length(y) - nrParams ;

4 % Estimate the standard deviation of y

5 stdY = norm(r) / sqrt(df);

6 % Evaluate the covariance

7 Bp = pinv( B );

8 Cov = stdY^2 * ( Bp * Bp’ ) ;

Listing 4: Code required to compute the covariances.

When considering the cases with other basis functions for the aperiodic
component, e.g. exponentials [14, 21] or polynomials [22], only line 5 in
Listing 1 and line 2 in Listing 3 need to be adjusted, according to the
appropriate basis functions. The rest of the code remains the same.

##  VIII Further extension

The approach described here could be further generalized to approximate
functions using rational polynomials, i.e. functions of the form

$r\left(x,\alpha,\beta\right)=\frac{f\left(x,\alpha\right)}{g\left(x,\beta\right)}.$
(29)

Such functions are commonly encountered in IIR-Filter design, in that case
polynomials in $z$. Zolotarev [32] found solutions to two specific cases based
on elliptical functions, but no general solutions. There are also solutions
based on Yule-Walker methods [33]. However, none of these previous solutions
take advantage of the method of variable projection to reduce the
dimensionality of the nonlinear optimization.

Let $\bm{B}_{f}(\bm{x},\bm{\alpha}),\bm{B}_{g}(\bm{x},\bm{\beta})$ be
polynomial basis matrices, such that the polynomials can be discretely written
as

$f\left(x,\bm{\alpha}\right)=\bm{B}_{f}(\bm{x},\bm{\alpha})\bm{\alpha},\quad
g\left(x,\bm{\beta}\right)=\bm{B}_{g}(\bm{x},\bm{\beta})\bm{\beta}.$ (30)

Then one could write

$r\left(x,\bm{\alpha},\bm{\beta}\right)=\underbrace{\left\\{\bm{B}_{g}(\bm{x},\bm{\beta})\bm{\beta}\right\\}^{-}\bm{B}_{f}(\bm{x},\bm{\alpha})}_{\bm{M}(\bm{x},\bm{\beta})}\bm{\alpha}=\bm{M}(\bm{x},\bm{\beta})\bm{\alpha}$
(31)

where $\left\\{\bm{B}_{g}\bm{\beta}\right\\}^{-}$ denotes some generalized
inverse, which needs to be further investigated. So if the function which
needed to be approximated was $d\left(x\right)$ and its shape was similar to
the shape of a rational function, one would write

$d\left(x\right)\approx\bm{M}(\bm{x},\bm{\beta})\bm{\alpha}$ (32)

which is a linear combination of nonlinear basis functions contained in the
matrix $\bm{M}(\bm{x},\bm{\beta})$, and thus the method of variable projection
could be used.

##  IX Conclusions

It has been shown that the new method permits the characterization of periodic
portions of a signal embedded in complex aperiodic backgrounds, cases where
the IEEE-standard 1057 is inappropriate. The use of $C^{n}$ continuous
B-Splines for the background is compatible with the piecewise local
interactions of a control system with an industrial process. It can be
concluded that the new method is a contribution to improving the
instrumentation of industrial processes. This new implementation works
especially well when the periodic portion is not perfectly periodic in its
observation time, since it avoids the Gibbs error and spectral leakage
associated with classical Fourier techniques. It also avoids the undesirable
effects of having to apply windowing to limit leakage. The algebraic
formulation, also yielded a simple approach to calculating the covariance of
the linear coefficients, i.e., all but one of the coefficients. Consequently,
the method provides a computation means of dealing with uncertainty. The
m-code snippets show the simplicity of the coding required to implement the
complete analysis.

## References

  * [1] P. O’Leary and D. Ninevski, “Estimating parameters of a sine wave by the method of variable projection,” in _2021 IEEE International Instrumentation and Measurement Technology Conference (I2MTC)_ , 2021, pp. 1–6.
  * [2] G. H. Golub and V. Pereyra, “The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate,” _SIAM Journal on Numerical Analysis_ , vol. 10, no. 2, pp. 413–432, 1973.
  * [3] G. Golub and V. Pereyra, “Separable nonlinear least squares: the variable projection method and its applications,” _Inverse Problems_ , vol. 19, pp. R1–R26(1), 01 2003.
  * [4] IEEE, “IEEE standard for digitizing waveform recorders,” _IEEE Std 1057-2017 (Revision of IEEE Std 1057-2007)_ , pp. 1–0, 2018.
  * [5] K. Hejn and A. Pacut, “Effective resolution of analog to digital converters,” _IEEE Instrumentation Measurement Magazine_ , vol. 6, no. 3, pp. 48–55, 2003.
  * [6] A. Baccigalupi, M. D’Arco, and A. Liccardo, “Parameters and methods for adcs testing compliant with the guide to the expression of uncertainty in measurements,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 66, no. 3, pp. 424–431, 2017.
  * [7] A. Albarbar, S. Mekid, A. Starr, and R. Pietruszkiewicz, “Suitability of mems accelerometers for condition monitoring: An experimental study,” _Sensors_ , vol. 8, no. 2, pp. 2192–2196, 2008.
  * [8] J. S. Lee, S. Choi, S. Kim, C. Park, and Y. G. Kim, “A mixed filtering approach for track condition monitoring using accelerometers on the axle box and bogie,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 61, no. 3, pp. 749–758, 2012.
  * [9] H. Wang, Z. Liu, A. Núñez, and R. Dollevoet, “Identification of the catenary structure wavelength using pantograph head acceleration measurements,” in _2017 IEEE International Instrumentation and Measurement Technology Conference (I2MTC)_ , 2017, pp. 1–6.
  * [10] M. Preeti, Koushik Guha, K. Baishnab, K. Dusarlapudi, and K. Narasimha Raju, “Low frequency mems accelerometers in health monitoring – a review based on material and design aspects,” _Materials Today: Proceedings_ , vol. 18, pp. 2152 – 2157, 2019, 2nd International Conference on Applied Sciences and Technology (ICAST-2019): Material Science. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S2214785319320310
  * [11] A. Jerri, _The Gibbs Phenomenon in Fourier Analysis, Splines and Wavelet Approximations_. Dordrecht, Netherlands: Kluwer Academic Publishers, 1998.
  * [12] F. Harris, “On the use of windows for harmonic analysis with the discrete fourier transform,” _Proceedings of the IEEE_ , vol. 66, pp. 55–83, 1978\.
  * [13] E. O. Brigham, _The Fast Fourier Transform and Its Applications_. USA: Prentice-Hall, Inc., 1988.
  * [14] A. R. Borden and B. C. Lesieutre, “Variable projection method for power system modal identification,” _IEEE Transactions on Power Systems_ , vol. 29, no. 6, pp. 2613–2620, 2014.
  * [15] G.-Y. Chen, M. Gan, C. L. P. Chen, and L. Chen, “A two-stage estimation algorithm based on variable projection method for gps positioning,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 67, no. 11, pp. 2518–2525, 2018.
  * [16] J. Abel, “A variable projection method for additive components with application to gps,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 30, no. 3, pp. 928–930, 1994.
  * [17] D. Robinson, S. Farsiu, and P. Milanfar, “Optimal Registration Of Aliased Images Using Variable Projection With Applications To Super-Resolution,” _The Computer Journal_ , vol. 52, no. 1, pp. 31–42, 04 2007. [Online]. Available: https://doi.org/10.1093/comjnl/bxm007
  * [18] J. Sheng and L. Ying, “A variable projection approach to parallel magnetic resonance imaging,” in _2008 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro_ , 2008, pp. 1027–1030.
  * [19] A. Aravkin, T. van Leeuwen, and N. Tu, “Sparse seismic imaging using variable projection,” in _ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2012.
  * [20] F. Gao, M. Li, A. Abubakar, and T. M. Habashy, “Application of variable projection scheme for data calibration in electromagnetic wave inversion,” in _2014 IEEE Antennas and Propagation Society International Symposium (APSURSI)_ , 2014, pp. 655–656.
  * [21] G.-Y. Chen, M. Gan, C. L. P. Chen, and H.-X. Li, “A regularized variable projection algorithm for separable nonlinear least-squares problems,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 2, pp. 526–537, 2019\.
  * [22] J. Handler, D. Ninevski, and P. O’Leary, “Decomposition of a periodic perturbed signal with unknown perturbation frequency by the method of variable projection,” in _7th International Conference on Mechanical Engineering and Automation Science (ICMEAS 2021)_ , (to appear) 2021.
  * [23] G. Wahba, _Spline Models for Observational Data_. Society for Industrial and Applied Mathematics, 1990. [Online]. Available: https://epubs.siam.org/doi/abs/10.1137/1.9781611970128
  * [24] L. Piegl and W. Tiller, _The NURBS Book_ , 2nd ed. New York, NY, USA: Springer-Verlag, 1996.
  * [25] J. Handler, September, 2021, MATLAB Central File Exchange, [Online]. https://de.mathworks.com/matlabcentral/profile/authors/23787014.
  * [26] G. Gowan and S. Brandt, _Data Analysis: Statistical and Computational Methods for Scientists and Engineers_ , ser. Ohlin Lectures; 7. Springer New York, 1998.
  * [27] K. Weierstrass, “über die analytische darstellbarkeit sogenannter willkürlicher functionen einer reellen veränderlichen,” _Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften zu Berlin, 1885 (II)_ , pp. 633–639, 789–805, 1885.
  * [28] C. de Boor, _A Practical Guide to Spline_ , 01 1978, vol. Volume 27.
  * [29] V. T. Dung and T. Tjahjowidodo, “A direct method to solve optimal knots of b-spline curves: An application for non-uniform b-spline curves fitting,” in _PloS one_ , 2017.
  * [30] D. Ninevski and P. O’Leary, “Detection of derivative discontinuities in observational data,” in _Advances in Intelligent Data Analysis XVIII_ , M. R. Berthold, A. Feelders, and G. Krempl, Eds. Cham: Springer International Publishing, 2020, pp. 366–378.
  * [31] D. O’Leary and B. Rust, “Variable projection for nonlinear least squares problems,” _Computational Optimization and Applications_ , vol. 54, no. 3, pp. 579–593, 2013.
  * [32] M.-P. Istace and J.-P. Thiran, “On the third and fourth zolotarev problems in the complex plane,” _SIAM Journal on Numerical Analysis_ , vol. 32, no. 1, pp. 249–259, 1995. [Online]. Available: http://www.jstor.org/stable/2158295
  * [33] B. Friedlander and P. Boaz, “The modified yule-walker method of arma spectral estimation,” _IEEE Transactions on Aerospace Electronic Systems_ , vol. AES-20, no. 2, 1984.


# Concentration of Hitting Times
in Erdős-Rényi graphs

Andrea Ottolini  Department of Mathematics, University of Washington, Seattle,
WA 98195, USA<EMAIL_ADDRESS>and  Stefan Steinerberger  Department of
Mathematics, University of Washington, Seattle, WA 98195, USA<EMAIL_ADDRESS>
###### Abstract.

We consider Erdős-Rényi graphs $G(n,p)$ for $0<p<1$ fixed and
$n\rightarrow\infty$ and study the expected number of steps, $H_{wv}$, that a
random walk started in $w$ needs to first arrive in $v$. A natural guess is
that an Erdős-Rényi random graph is so homogeneous that it does not really
distinguish between vertices and $H_{wv}=(1+o(1))n$. Löwe-Terveer established
a CLT for the Mean Starting Hitting Time suggesting
$H_{wv}=n\pm\mathcal{O}(\sqrt{n})$. We prove the existence of a strong
concentration phenomenon: $H_{wv}$ is given, up to a very small error of size
$\lesssim\sqrt{\log{n}}/\sqrt{n}$, by an explicit simple formula involving
only the total number of edges $|E|$, the degree $\deg(v)$ and the distance
$d(v,w)$.

###### Key words and phrases:

Erdős-Rényi graphs, hitting time, random walk.

######  Mathematics Subject Classification:

60J10, 05C80, 60B20

S.S. was supported by the NSF (DMS-2123224) and the Alfred P. Sloan
Foundation.

##  1\. Introduction and Results

###  1.1. The Phenomenon.

The purpose of this paper is to demonstrate a very strong concentration
phenomenon for hitting times on Erdős-Rényi random graphs. Given a realization
of the graph, the hitting time $H_{wv}$ is the expected time that a simple
random walk started in vertex $w$ needs to reach the vertex $v$ for the first
time. We consider the distribution of hitting times $H_{wv}$ when the graph
$G=G(n,p)$ is an Erdős-Rényi random graph for $0<p<1$ fixed and
$n\rightarrow\infty$. The phenomenon we are interested in can be illustrated
with a simple example, see Fig. 1.

Figure 1. Left: the distribution of all hitting times $H_{ij}$ for a
realization of $G(1000,1/2)$. Right: the distribution of hitting times
$H_{wv}$ for $w\neq v$ for a fixed vertex $v$. This quantity is tightly
concentrated in two very short intervals.

A natural first guess is that such a random walk is rapidly mixing and should
spend roughly equal amounts of time in each vertex and all the vertices
should, at least to leading order, be indistinguishable. This would suggest
that $H_{wv}\sim n$ whenever $w\neq v$ and this was also suggested in the
physics literature [13]. A version of this scaling result was rigorously
proven by Löwe & Torres [9]. If we denote by $H_{\mu v}:=\sum\mu(w)H_{wv}$ the
hitting time starting from a point distributed according to a probability
measure $\mu$, then $H_{\pi v}=(1+o(1))n$ where $\pi(w):=\deg(w)/(2E)$ is the
(random) stationary distribution. We also refer to subsequent results of
Helali & Löwe [4], von Luxburg, Radl & Hein [11] and Sylvester [14]. A look at
some numerical examples (see Fig. 1) suggests that the typical deviation of
$H_{wv}$ from the mean appears to be on the order of $\sqrt{n}$. One would
naturally expect the existence of a central limit theorem. Such a result was
recently established by Löwe & Terveer [10], who showed that $H_{\pi v}-n$,
suitably rescaled, converges to a Gaussian (in distribution). This, while not
a statement about $H_{wv}$ in itself, does indicate that averaged hitting
times have Gaussian fluctuations.

###  1.2. Main Result.

We provide what is essentially an explicit formula for $H_{wv}$ up to a very
small error term that tends to 0 as $n\rightarrow\infty$.

######  Theorem.

Let $G(n,p)$ be an Erdős-Rényi random graph with $0<p<1$. Then, as
$n\rightarrow\infty$, we have with high probability that for any vertex $v$
and all vertices $w\neq v$

$H_{wv}=\frac{2|E|}{\deg(v)}+\begin{cases}-1\quad&\mbox{if}~{}(w,v)\in E\\\
-1+1/p\quad&\mbox{if}~{}(w,v)\notin
E\end{cases}\quad+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)$

where the implicit constant in the error term depends only on $p$.

This explains the second picture in Fig. 1, we refer to Fig. 2 for more
examples. We also prove (Lemma 1) that the average hitting time starting from
vertices $w$ that are adjacent to $v$ is exactly $2|E|/\deg(v)-1$. The Theorem
implies that the central limit theorem obtained by Löwe-Terveer [10] for
$H_{\pi\cdot}$, the hitting time from a stationary initial point, holds for
$H_{\mu v}$ for any starting measure $\mu$, and in particular it also holds
for the hitting times $H_{wv}$ themselves, see §2.8.

Figure 2. Left: $G(4000,0.2)$ with predicted localization at
$2|E|/\deg(v)-1=3965.3$ and $2|E|/\deg(v)+4=3970.3$. Right: $G(4000,0.8)$ with
predicted localization at $3994.9$ and $3996.16$.

One interpretation, which also informs the structure of the proof, is that an
Erdős-Rényi random graph has diameter 2 (with high probability) and can be
understood as behaving a little bit like a three-state Markov chain: the
vertex $v$ under consideration, the vertices at distance 1 from $v$, and the
vertices at distance 2 from $v$. Note that transitions to $v$ are rare, given
that a typical vertex has an high degree. During that time, the mixing
properties of the graph are of sufficiently high quality that they essentially
lead to near-deterministic behavior. However, it is important whether one
starts at distance 1 or 2 from the target vertex $v$ (not surprising since we
aim to determine the hitting term up to an error of $o(1)$). In fact, since
transitions from most vertices at distance $2$ to the set of vertices at
distance $1$ happen at rate roughly $\sim p$, the correction $1/p$ can be
understood as a geometric hitting time. It is conceivable that the result has
an extension when $p$ decays with $n$, generalizing accordingly the number of
possible states since the diameter will increase with high probability (see,
e.g., [1] for a variety of regimes). Moreover, while some features of our
proof rely on the Erdős-Rényi structure, a similar approach may lead to less
refined results (e.g., a law of large number or a central limit theorem) for
other models of random graphs, provided that the mixing time is sufficiently
small compared to the hitting time, bypassing a spectral approach. This seems
like a promising avenue for further research.

###  1.3. Implications.

The result has a number of implications in terms of known results. Löwe &
Torres [9] prove, using a spectral theory approach, that the averaged hitting
time $H_{\pi v}$ satisfies

$H_{\pi v}\geq\frac{2|E|}{\deg(v)}-2$

which can now be seen to be almost optimal. Our Theorem implies

$\displaystyle H_{\pi
v}=\sum_{w}\pi(w)H_{wv}=\frac{2|E|}{\deg(v)}-2+\frac{1}{p}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right).$

The result also has immediate implications for the commute time
$\kappa_{wv}=H_{wv}+H_{vw}$ and the effective resistance

$\sigma_{wv}=\frac{\kappa_{wv}}{2|E|}.$

For example, we can deduce a very precise asymptotic expansion

$\sigma_{wv}=\frac{1}{\deg(v)}+\frac{1}{\deg(w)}+\frac{2}{n^{2}p}\begin{cases}-1\quad&\mbox{if}~{}(v,w)\in
E\\\ -1+1/p\quad&\mbox{if}~{}(v,w)\notin
E\end{cases}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{n^{5/2}}\right).$

We remark that the effective resistance is intimately connected to a number of
classical objects in probability, such as uniform spanning trees (a result
dating back to Kirchhoff [6]), and it is easy to rephrase the asymptotic above
correspondingly (see, e.g., Chapter $4$ in [12]). Comparing to arbitrary
graphs, a result of Lovász [8] shows that the commute time always satisfies

$|E|\left(\frac{1}{\deg(v)}+\frac{1}{\deg(w)}\right)\leq\kappa_{wv}\leq\frac{2|E|}{1-\lambda_{2}}\left(\frac{1}{\deg(v)}+\frac{1}{\deg(w)}\right),$

where $\lambda_{2}$ is the second eigenvalue of the transition matrix. For a
typical Erdős-Rényi random graph, our result says that the lower bound is
almost saturated. We conclude with some spectral implication (see Löwe &
Torres [9] or Lovász [8] for additional details). Let
$A\in\left\\{0,1\right\\}^{n\times n}$ denote the adjacency matrix of the
graph and $D$ denote the degree matrix of the graph, we can introduce
$B=D^{-1/2}AD^{-1/2}$. $B$ has largest eigenvalue $\lambda_{1}=1$ and has
eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{n}$. We use
$\sigma_{k}$ to denote the eigenvectors of length 1 associate to the $k-$th
eigenvalue. Then, see Löwe & Torres [9] or Lovász [8], we have

$H_{\pi
v}=\frac{2|E|}{\deg(v)}\sum_{w=2}^{n}\frac{\sigma_{wv}^{2}}{1-\lambda_{k}}.$

Since $H_{\pi v}$ is nearly constant in $v$, this can be understood as a type
of weighted equidistribution result for the eigenvectors of $B$. Finally, our
result allows for estimates related to the quasi-stationary distribution
$\tilde{\pi}_{v}$ associated to a vertex $v$, that is to say, the Perron-
Frobenius eigenvector of the sub-stochastic matrix obtained by neglecting the
row and column corresponding to $v$ in the transition matrix (see [3] and
references therein for a discussion on its use). It is a classical fact that
its eigenvalue $\lambda_{v}$ satisfies

$\displaystyle H_{\tilde{\pi}_{v}v}=\frac{1}{1-\lambda_{v}},$

and so our result gives two-sided bounds on $\lambda_{v}$ that are uniform in
$v\in V$ (w.h.p.).

##  2\. Proof

###  2.1. Outline

It is a classical result (see, for example, [7]) that with high probability
the diameter of $G(n,p)$ for fixed $p$ is 2. From now on, we will work on the
event that this occurs. We fix a vertex $v\in V$ and consider from now on only
the (first) hitting time of $v$ of a random walk starting in any vertex $w\in
V$. Naturally, $H_{vv}=0$. For any other vertex $w\neq v$, we have

$\displaystyle H_{wv}=1+\frac{1}{\deg(w)}\sum_{(w,u)\in E}H_{uv}.$ (1)

Having fixed $v\in V$, we denote the set of all vertices at distance 1 from
$v$ by $A$,

$A=\left\\{w\in V:d(v,w)=1\right\\}$

and the set of all vertices at distance 2 from $v$ by
$B=V\setminus\left(\left\\{v\right\\}\cup A\right)$. We first establish
concentration of the hitting times in $A$ and $B$ and then use the
concentration to deduce everything else. The proof works as follows.

  1. (1)

We first recall some basic facts from probability theory (in §2.2) that will
be used throughout to argument: concentration of the degrees, results on the
spectral gap, and an exact computation for the hitting time when the starting
distribution is uniform in $A$.

  2. (2)

In §2.3 we derive a combinatorial formula for the hitting time. It is then
shown that many of the terms in the formula are tightly concentrated,
conditional on a certain a priori bound to be established.

  3. (3)

§2.4 shows the a priori bound

$H_{wv}=n+\mathcal{O}(\sqrt{n\log{n}})$

with high probability. This result is not as good as the main result but shows
that certain terms that arise in the formula from §2.3 are more tightly
concentrated than one might expect from $H_{wv}=(1+o(1))\cdot n$.

  4. (4)

§2.5 uses all the preceeding arguments to show strong concentration of the
hitting time in the set $B$,

$\max_{w_{1},w_{2}\in
B}|H_{w_{1}v}-H_{w_{2}v}|\lesssim\frac{\sqrt{\log{n}}}{\sqrt{n}}.$

  5. (5)

In §2.6 a similar argument is repeated for vertices in $A$ and

$\max_{w_{1},w_{2}\in
A}|H_{w_{1}v}-H_{w_{2}v}|\lesssim\frac{\sqrt{\log{n}}}{\sqrt{n}}.$

  6. (6)

Having shown concentration in both $A$ and $B$, it is shown in §2.7 that the
difference in expectation between the hitting times in $A$ and $B$ is
$1/p+\mathcal{O}(\sqrt{\log{n}}/\sqrt{n})$. Since we already know (Lemma 1)
the expected hitting time for vertices in $A$, the result follows.

###  2.2. Preliminary facts

We first collect a couple of elementary facts that will be used several times.
We will often omit w.h.p. (with high probability) in what follows: every time
we invoke a property of Erdős-Rényi random graphs it is to be understood as
being valid with high probability. The sequence of degrees are Bernoulli
random variables $\deg(v)\sim B(n,p)$, with $\mathbb{E}\deg(v)=np$. While they
are not independent, a union bound, together with a Chernoff bound for
binomial random variables, gives

$np-c_{p}\sqrt{n}\sqrt{\log{n}}\leq\min_{v\in V}\deg(v)\leq\max_{v\in
V}\deg(v)\leq np+c_{p}\sqrt{n}\sqrt{\log{n}}$

We note that more precise results, in particular about the size of $c_{p}$,
are known (see, e.g., [2]), though we will not try to optimize over that as
they are not required for the remainder of the argument. In particular, this
gives

$\max_{v\in V}\left|\deg(v)-np\right|\lesssim\sqrt{n\log{n}}$ (2)

with high likelihood and with an implicit constant depending on $p$ only.
Another ingredient that we will use several times is that the spectral gap of
the random walk on an Erdős-Rényi random graph is, with high probability is
$\leq c_{p}n^{-1/2}$. Once more, we will not try to optimize over the
dependence on $p$, though the interested reader should refer to [5]. In
particular, this implies that, with high probability, one has a bound in total
variation between $\nu_{k}$, the distribution of a walker with initial
distribution $\nu$ after $k$ steps, and its stationary distribution $\pi$. The
Cauchy-Schwarz inequality together with an upper bound on the spectral gap
shows that

$\displaystyle\|\nu_{k}-\pi\|_{1}$
$\displaystyle\leq\sqrt{n}\cdot\|\nu_{k}-\pi\|_{2}$ (3)
$\displaystyle\lesssim\sqrt{n}\cdot|\lambda_{2}|^{k}\cdot\|\nu_{0}-\pi\|_{2}\lesssim
c_{p}^{k}\cdot n^{\frac{1-k}{2}}.$

This is consistent with the fact that after one step the total variation
distance between $\nu_{1}$ and $\pi$ bounded from below if $\nu_{0}$ is a
point mass – $\nu_{1}$ is supported on approximately $pn$ vertices – while
after two steps mixing already occurs. Since hitting times have a size
$\mathcal{O}(n)$, it becomes clear that considering three steps of the Markov
chain is typically needed to get results that are precise up to constant
order.

Our argument will benefit from having at least one probability measure $\nu$
for which we can obtain explicit estimates on the expected hitting time.

######  Lemma 1.

Consider the probability measure $\nu=1_{A}/\deg(v)$. Then, we have

$\displaystyle H_{\nu v}=\frac{2|E|}{\deg(v)}-1.$

###### Proof.

Consider a random walker starting from the stationary distribution $\pi$ given
by $\pi(w)=\deg(w)/(2|E|)$. There are two cases: either the initial vertex is
already $v$ (and then the hitting time is 0) or it is not, in which case the
new probability distribution after one step is distributed according to some
new measure $\mu$. On the other hand, conditional on the starting point being
$v$, the distribution after one step would be exactly $\nu$. Therefore, we
have

$\displaystyle\mu=\frac{\pi-\pi(v)\nu}{1-\pi(v)}.$

In particular, we obtain

$\displaystyle H_{\pi v}=(1-\pi(v))\left(1+H_{\mu v}\right)=1-\pi(v)+H_{\pi
v}-\pi(v)H_{\nu v}.$

Rearranging, this gives

$\displaystyle H_{\nu v}=\frac{1}{\pi(v)}-1=\frac{2|E|}{\deg(v)}-1.$

∎

The result has an electrical network interpretation, which can be turned into
a proof when $A=\\{w\\}$ is a single vertex connected to $v$ through $\deg(v)$
edges. The effective resistance $\kappa_{wv}$ is that of $\deg(v)$ resistances
in parallel. Then,

$\displaystyle\frac{H_{wv}+H_{vw}}{2E}=\kappa_{vw}=\frac{1}{\deg(v)},$

and the result follows since $H_{vw}=1$.

###  2.3. A formula for the Hitting time for $w\in B$.

Let us now fix a vertex $w\in B$. As an initial motivation, we consider the
distribution of a random walk started in $w\in B$ with a length of 2 steps. We
start with the basic fact (introduced above as (1)) that for a vertex $w\neq
v$

$\displaystyle H_{wv}=1+\frac{1}{\deg(w)}\sum_{(w,u)\in E}H_{uv}.$

Since $w\in B$ and $d(w,v)=2$, we can apply the identity twice. If we denote
the arising distribution after 2 steps by $\nu_{w}$, then this implies the
exact equation

$H_{wv}=2+H_{\nu_{w}v}.$

For the remainder of the argument, we will consider the distribution of random
walks started in $w\in B$ after 3 steps and this probability distribution will
be denoted by $\mu_{w}$. The analogous equation

$\displaystyle H_{wv}=3+H_{\mu_{w}v}\qquad\mbox{is false}$ (4)

because (1) is only valid in vertices $w\neq v$ but a random walk started in
$w$ might end up in $v$ after 2 steps. We also note that this is the only
source of error. If the random walk ends up in $v$ after 2 steps, then it is
uniformly distributed in $A$ after 3 steps and the equation (4) mistakenly
contributes

$\mathbb{P}\left(\mbox{random walk is in}~{}v~{}\mbox{after 2
steps}\right)\cdot\frac{1}{|A|}\sum_{w\in A}H_{wv}$

while the true hitting time would have been 2. We note that the sum, the
average hitting time for a starting vertex in $A$, has already been computed
above in Lemma 1 but we will not actually need its precise value at this point
in the argument. Correcting for this, we arrive at

$\displaystyle H_{wv}$ $\displaystyle=3+\sum_{a\in V}\mu_{w}(a)\cdot H_{av}$
(5) $\displaystyle+\mathbb{P}\left(\mbox{random walk is in}~{}v~{}\mbox{after
2 steps}\right)\cdot\left(-1-\frac{1}{|A|}\sum_{a\in A}H_{av}\right).$

We will now argue that (5) can be used to show that the hitting time is
tightly concentrated by showing that each of its terms is tightly
concentrated. We start by noting that Lemma 1 implies that

$\left(-1-\frac{1}{|A|}\sum_{a\in A}H_{av}\right)=-\frac{2|E|}{\deg(v)}$

is independent of $w$. We continue by computing the probability of a random
walk started in $w$ hitting $v$ within 2 steps. For this to happen, the first
step of the random walk has to lead to $A$ and the second step has to lead
from $A$ to $v$. The likelihood of moving from $w$ to $A$ is simply the number
of neighbors that $w$ has in $A$ compared to its total degree. Introducing the
abbreviation

$N_{A}(w)=\left\\{a\in A:(a,w)\in E\right\\}$

for the number of neighbors a vertex has in $A$, we have that the likelihood
of moving from $w$ to $A$ is given by

$\frac{\\#N_{A}(w)}{\deg(w)}=\frac{|A|}{n}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)\qquad\mbox{with
high probability.}$

owing to (2). The second ingredient is the likelihood of moving from $A$ to
$v$ in the second step which, conditioned on the first step leading from $w$
to $A$, happens with likelihood

$\frac{1}{\\#N_{A}(w)}\sum_{z\in N_{A}(w)}\frac{1}{\deg(z)}.$

This is an average over inverse degrees over the set of neighbors of $w$ that
are in $A$. Abbreviating

$X=\left|\frac{1}{\\#N_{A}(w)}\sum_{z\in
N_{A}(w)}\frac{1}{\deg(z)}-\frac{1}{(n-1)p}\right|$

we obtain, appealing (2) that (w.h.p.)

$\displaystyle X$ $\displaystyle\leq\max_{v\in
V}\left|\frac{1}{\deg(v)}-\frac{1}{(n-1)p}\right|$
$\displaystyle\leq\left|\frac{1}{np\pm\mathcal{O}(\sqrt{n\log{n}})}-\frac{1}{np}\right|\lesssim\frac{\sqrt{n\log{n}}}{n^{2}p^{2}}\lesssim\frac{\sqrt{\log{n}}}{n^{3/2}}.$

We see that $\mathbb{P}=\mathbb{P}\left(\mbox{random walk
in}~{}v~{}\mbox{after 2 steps}\right)$ satisfies

$\displaystyle\mathbb{P}$
$\displaystyle=\left(\frac{|A|}{n}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)\right)\left(\frac{1}{np}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{n^{3/2}}\right)\right)$
(6)
$\displaystyle=\frac{|A|}{n^{2}p}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{n^{3/2}}\right).$

###  2.4. A Cheap Uniform Bound

The purpose of this short section is to deduce a cheap uniform bound. This
estimate will then be bootstrapped in the subsequent sections to obtain the
main result.

######  Lemma 2.

Let $G=G(n,p)$ with $0<p<1$. Then, as $n\rightarrow\infty$, with high
probability

$\max_{v\neq w\in V}~{}\left|H_{wv}-n\right|\lesssim\sqrt{n\log{n}}.$

and also

$\max_{v\in V,w_{1},w_{2}\in V\setminus\\{v\\}}|H_{w_{1}v}-H_{w_{2}v}|\lesssim
1.$

In both results, the implicit constants depend on $p$ only.

###### Proof.

We know from (2) that, with high likelihood,

$\max_{w\in V}\left|\deg(w)-np\right|\lesssim\sqrt{n\log{n}}.$

We first argue that $H_{wv}\leq C_{p}n$ for some constant $0<C_{p}<\infty$.
Let $X_{0},X_{1},\ldots,$ denote the simple random walk on the graph. We note
that there are $|A|=np\pm\mathcal{O}(\sqrt{n\log{n}})$ vertices where the
likelihood of hitting $v$ in the next step is at least

$\mathbb{P}\left(X_{k+1}=v\big{|}X_{k}\in
A\right)=\frac{1}{\deg(X_{k})}=\frac{1}{np}\pm\mathcal{O}\left(\frac{\sqrt{\log{n}}}{n^{3/2}}\right).$

If the random walk $X_{k}\in B$, then

$\mathbb{P}\left(X_{k+1}\in A\big{|}X_{k}\in B\right)\geq p+o(1).$

We can now do a simple case distinction: either $X_{k}\in A$ in which case
$X_{k+1}=v$ with likelihood at least $1/(np)$ or $x_{k}\in B$ in which case
$X_{k+2}=v$ with likelihood at least $1/n$. Using a geometric probability
distribution for stochastic domination, we see that $H_{wv}\leq C_{p}n$. We
can now consider an arbitrary starting vertex $w\in V$ together with a random
walk of 3 steps started in $w$. Let $T_{wv}$ be the time that it takes to hit
$v$ starting from $w$, so that $\mathbb{E}[T_{wv}]=H_{wv}$. There are two
cases: the random walk happens to hit $v$ within the first $3$ steps or it
does not. Therefore, if $\mu_{3}$ denote the distribution after $3$ steps,

$\displaystyle H_{wv}$ $\displaystyle=\mathcal{O}(1)+\sum_{u\in
V}\mu_{3}(u)\mathbb{E}[T_{uv}\mathbb{E}[1_{T_{wv}>3}|X_{3}=u]].$

Note that

$\mathbb{E}[T_{uv}\mathbb{E}[1_{T_{wv}>3}|X_{3}=u]]\leq H_{uv}\leq C_{p}n.$

The same bound holds if we replace $w$ with the stationary distribution. Using

$\displaystyle\|\mu_{3}-\pi\|_{1}\lesssim\frac{1}{n},$

which follows from (3) with $k=3$ we deduce that with high probability and
uniformly over all $v\in V$ and $w_{1},w_{2}\in V\setminus\\{v\\}$

$|H_{w_{1}v}-H_{w_{2}v}|\lesssim 1.$

Note that, by convexity, the result remains true if $w_{2}$ is replaced by an
arbitrary probability measure. If we choose $\nu$ from Lemma 1, we deduce

$\displaystyle\left|H_{wv}-\frac{2E}{\deg j}+1\right|\lesssim 1,$

and thus we deduce since, using (2), that

$\displaystyle\left|\frac{2E}{\deg j}-n\right|=\mathcal{O}\left(\sqrt{n\log
n}\right)$

with high probability. ∎

###  2.5. Concentration of Hitting Times in $B$.

We can now use these results to establish concentration of the hitting times
in $B$. For this purpose, let us keep our original arbitrary vertex $w\in B$
considered above and let us additionally consider another vertex $w_{2}\in B$.
The distribution of a random walk started in $w_{2}$ after 3 steps will be
denoted by $\mu_{2}$. Using (5) once in $w$ and once in $w_{2}$ and invoking
the uniform asymptotics (6) we arrive at

$\displaystyle H_{wv}-H_{w_{2}v}$
$\displaystyle=\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+\sum_{c\in
V}(\mu(c)-\mu_{2}(c))\cdot H_{cv}$

Lemma 2, the fact that two probability distributions $\mu$ and $\mu_{2}$ have
equal $\ell^{1}-$norm and the Cauchy-Schwarz inequality allow us to rewrite
this as

$\displaystyle\left|H_{wv}-H_{w_{2}v}\right|$
$\displaystyle=\left|\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+\sum_{c\in
V}(\mu(c)-\mu_{2}(c))\cdot(n+(H_{cv}-n))\right|$
$\displaystyle\leq\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+\sum_{c\in
V}\left|\mu(c)-\mu_{2}(c)\right|\cdot\left|H_{cv}-n\right|$
$\displaystyle\lesssim\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+\sqrt{n\log{n}}\cdot\|\mu-\mu_{2}\|_{\ell^{1}}.$

Using a triangular inequality and (3) with $k=3$, we obtain

$\|\mu-\mu_{2}\|_{\ell^{1}}\lesssim\frac{1}{n},$

leading to

$\displaystyle\left|H_{wv}-H_{w_{2}v}\right|$
$\displaystyle\lesssim\frac{\sqrt{\log{n}}}{\sqrt{n}}.$ (7)

###  2.6. Concentration for vertices in $A$.

We will now adapt the argument for vertices in $A$. Let $w\in A$ and consider
the usual random walk in $A$ after 3 steps. We would like to argue analogously
as before to obtain a combinatorial formula for the hitting time. There are
now two types of mistakes that can happen: those that arise if the random walk
already ends up in $v$ after 1 step or after 2 steps.

  1. (1)

Scenario 1. The random walk moves from $w$ to $v$ in the first step. In that
case the true hitting time is 1 where as the naive formula (4) produces an
error given by the average hitting time in a point after 2 random steps
started in $v$ (as before, this quantity is $\mathcal{O}(n)$ and independent
of $w$).

  2. (2)

Scenario 2. The random walk arrives in $v$ at the second step (but not the
first). In that case, we deduce that the first random step has to lead from
$w$ to another vertex in $A$ and then from that vertex to $v$. The error made
in the naive formula is that it produces the average hitting time in $A$ as
opposed to the true hitting time, which equals to $2$.

We deduce that there are two numbers $O_{1},O_{2}=\mathcal{O}(n)$, independent
of $w\in A$, such that

$\displaystyle H_{wv}$ $\displaystyle=3+\sum_{a\in V}\mu(a)\cdot H_{av}$ (8)
$\displaystyle+\mathbb{P}\left(\mbox{Scenario 1}\right)\cdot
O_{1}+\mathbb{P}\left(\mbox{Scenario 2}\right)\cdot O_{2}.$

Note that $O_{1}$ is simply the expected hitting time weighted by the outcome
of a random walk of length 2 started in $v$ and $O_{2}$ is the expected
hitting time weighted by the distribution of a random walk of length 1 started
in $v$ (which, as above, is merely the average hitting time in $A$). It
remains to show that, just as above, the likelihoods of the two scenarios
depend on $w$ in a weak sense. We start with Scenario 1. We have, using (2),

$\mathbb{P}\left(\mbox{Scenario
1}\right)=\frac{1}{\deg(w)}=\frac{1}{np\pm\mathcal{O}(\sqrt{n\log{n}})}=\frac{1}{np}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{n^{3/2}}\right).$

The likelihood of Scenario 2 was already implicitly computed above: there we
computed the likelihood of having a random walk from a vertex in $B$ travel to
a vertex in $A$ and then to $v$. However, that likelihood is independent of
whether one starts in $A$ or in $B$ and the very same argument as above (see
(6)) gives

$\displaystyle\mathbb{P}\left(\mbox{Scenario
2}\right)=\frac{|A|}{n^{2}p}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{n^{3/2}}\right)$

The remainder of the argument is exactly the same as above and we deduce that

$\displaystyle\max_{w_{1},w_{2}\in
A}|H_{w_{1}v}-H_{w_{2}v}|\lesssim\frac{\sqrt{\log{n}}}{\sqrt{n}}.$ (9)

###  2.7. Proof of the Theorem

At this point, we know that the expected hitting time in both $A$ and $B$ is
essentially constant, owing to (9) and (7), up to an error of size
$\lesssim\sqrt{\log n}/{n}$. We define the average hitting times in $A$ and in
$B$

$\displaystyle H_{A}=\frac{1}{|A|}\sum_{w\in A}H_{wv},\quad
H_{B}=\frac{1}{|B|}\sum_{w\in B}H_{wv}.$

From Lemma (1), we know

$H_{A}=\frac{2E}{\deg j}-1,$ (10)

and, appealing to (9), we obtain that for $w\in A$

$\displaystyle H_{wv}=\frac{2E}{\deg j}-1+\mathcal{O}\left(\frac{\sqrt{\log
n}}{\sqrt{n}}\right).$

It remains to show that

$H_{B}-H_{A}=\frac{1}{p}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right),$

for the result will follow from (7) and (10). Let now $w\in B$ be an arbitrary
vertex. Performing one step of a random walk in $w$, we deduce that

$\displaystyle H_{B}$
$\displaystyle=H_{wv}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)$
$\displaystyle=\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+1+\frac{1}{\deg(w)}\sum_{z\in
N_{A}(w)}H_{zv}+\frac{1}{\deg(w)}\sum_{z\in N_{B}(w)}H_{zv}.$

At this point we can invoke the strong concentration of hitting times in $A$
and $B$ once more to deduce that

$\displaystyle H_{B}$
$\displaystyle=\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+1+\frac{1}{\deg(w)}\sum_{z\in
N_{A}(w)}H_{A}+\frac{1}{\deg(w)}\sum_{z\in N_{B}(w)}H_{B}$
$\displaystyle=\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+1+\frac{\\#N_{A}(w)}{\deg(w)}H_{A}+\frac{\\#N_{B}(w)}{\deg(w)}H_{B}.$

Vertices in $B$ only have neighbors in $A$ and in $B$ and therefore, for $w\in
B$,

$1=\frac{\\#N_{A}(w)+\\#N_{B}(w)}{\deg(w)}$

from which we deduce

$\frac{\\#N_{A}(w)}{\deg(w)}H_{B}=\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+1+\frac{\\#N_{A}(w)}{\deg(w)}H_{A}$

from which we deduce

$H_{B}=\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+\frac{\deg(w)}{\\#N_{A}(w)}+H_{A}.$

We also note that, from (2),

$\displaystyle\frac{\deg(w)}{\\#N_{A}(w)}$
$\displaystyle=\frac{np+\mathcal{O}(\sqrt{n\log{n}})}{p|A|+\mathcal{O}(\sqrt{n\log{n}})}$
$\displaystyle=\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right)+\frac{n}{|A|}=\frac{1}{p}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right).$

Thus

$\displaystyle
H_{B}-H_{A}=\frac{1}{p}+\mathcal{O}\left(\frac{\sqrt{\log{n}}}{\sqrt{n}}\right).$

###  2.8. A Central Limit Theorem

We conclude by showing how the Theorem immediately implies a central limit
theorem for the hitting times starting from any point. It can be easily
generalized to an arbitrary initial configuration that does not include the
target point. This follows from the Theorem and the fact that the main
fluctuation comes from the fluctuation of the degree of the target vertex
while being virtually independent (up to fluctuations of size
$\mathcal{O}(1)$) of everything else.

######  Corollary.

Fix two vertices $1,2$. Then, the distribution of $H_{12}$ over all Erdős-
Rényi random graphs satisfies, as $n\rightarrow\infty$,

$\displaystyle\sqrt{\frac{p}{n(1-p)}}\left(H_{12}-n\right)\Rightarrow\mathcal{N}(0,1)$

###### Proof.

The central limit theorem for a binomial random variable says that

$\displaystyle\deg(2)=np\left(1-\sqrt{\frac{1-p}{pn}}Z_{n}\right),$

where $Z_{n}$ converges in distribution to a normal random variable. A
standard concentration bound for the number of edges gives

$\displaystyle 2|E|=n^{2}p\left(1+\mathcal{O}\left(\frac{\sqrt{\ln
n}}{n}\right)\right)$

with high probability. Combining with our main theorem, we obtain

$H_{12}=n\frac{1+\mathcal{O}\left(\frac{\sqrt{\ln
n}}{n}\right)}{1-\sqrt{\frac{1-p}{np}}Z_{n}}+\mathcal{O}\left(1\right),$

where the error terms are uniform with respect to the measures $\mu_{n}$.
Therefore

$\displaystyle\sqrt{\frac{p}{n(1-p)}}\left(H_{12}-n\right)$
$\displaystyle=\sqrt{\frac{np}{1-p}}\left[\frac{1+\mathcal{O}\left(\frac{\sqrt{\ln
n}}{n}\right)}{1-\sqrt{\frac{1-p}{np}}Z_{n}}-1+\mathcal{O}\left(\frac{1}{n}\right)\right]$
$\displaystyle=\frac{Z_{n}+\mathcal{O}\left(\frac{\sqrt{\ln
n}}{\sqrt{n}}\right)}{1+\mathcal{O}\left(\frac{1}{\sqrt{n}}\right)}+\mathcal{O}(1).$

The result then follows from Slutsky’s theorem. ∎

Acknowledgment. The authors are grateful for discussions with Karel Devriendt.

## References

  * [1] F. Chung, L. Linyuan, The diameter of sparse random graphs. Advances in Applied Mathematics 26.4 (2001): 257-279.
  * [2] A. Desolneux, L. Moisan, and J. Morel. From gestalt theory to image analysis: a probabilistic approach. Vol. 34. Springer Science & Business Media, 2007.
  * [3] P. Diaconis, L. Miclo. On quantitative convergence to quasi-stationarity. Annales de la Faculte des sciences de Toulouse: Mathematiques. Vol. 24. No. 4. 2015.
  * [4] A. Helali and M. Löwe, Hitting times, commute times, and cover times for random walks on random hypergraphs. Statist. Probab. Lett. 154 (2019), 108535, 6 pp.
  * [5] C. Hoffman, M. Kahle, E. Paquette. Spectral gaps of random graphs and applications. International Mathematics Research Notices 2021.11 (2021): 8353-8404.
  * [6] G. Kirchhoff, Über die Auflösung der Gleichungen, auf welche man bei der Untersuchung der linearen Vertheilung galvanischer Ströme geführt wird. Ann. Phys. Chem. 72 (1847), 497–508
  * [7] V. Klee and D. Larman, Diameters of random graphs, Canadian Journal of Mathematics 33.3 (1981): 618-640.
  * [8] L. Lovász, Random walks on graphs: a survey. (English summary) Combinatorics, Paul Erdős is eighty, Vol. 2 (Keszthely, 1993), 353–397, Bolyai Soc. Math. Stud., 2, János Bolyai Math. Soc., Budapest, 1996.
  * [9] M. Löwe and F. Torres, On hitting times for a simple random walk on dense Erdős-Rényi random graphs. Statistics & Probability Letters 89 (2014), p. 81–88.
  * [10] M. Löwe and S. Terveer, A Central Limit Theorem for the average target hitting time for a random walk on a random graph. arXiv preprint arXiv:2104.01053.
  * [11] U. von Luxburg, A. Radl and M. Hein, Hitting and commute times in large random neighborhood graphs. The Journal of Machine Learning Research, 15 (2014), 1751–1798.
  * [12] R. Lyons, Y. Peres, Probability on trees and networks. Cambridge University Press, 2017.
  * [13] V. Sood, S. Redner and D. ben-Avraham, First-passage properties of the Erdős-Renyi random graph. J. Phys. A 38 (2005), no. 1, 109–123.
  * [14] J. Sylvester, John Random walk hitting times and effective resistance in sparsely connected Erdős-Rényi random graphs. J. Graph Theory 96 (2021), no. 1, 44–84.


# Kinematic Evidence of an Embedded Protoplanet in HD 142666 Identified by
Machine Learning

J. P. Terry Department of Physics and Astronomy, The University of Georgia,
Athens, GA 30602, USA.  Center for Simulational Physics, The University of
Georgia, Athens, GA 30602, USA.  C. Hall Department of Physics and Astronomy,
The University of Georgia, Athens, GA 30602, USA.  Center for Simulational
Physics, The University of Georgia, Athens, GA 30602, USA.  S. Abreau
Department of Medicine, Division of Cardiology, University of California, San
Francisco, CA, 94143, USA.  Cardiovascular Research Institute, San Francisco,
CA, 94158, USA.  S. Gleyzer Department of Physics and Astronomy, The
University of Alabama, Tuscaloosa, AL 35487, USA

(Received January 12, 2023; Revised February 17, 2023; Accepted March 10,
2023)

###### Abstract

Observations of protoplanetary disks have shown that forming exoplanets leave
characteristic imprints on the gas and dust of the disk. In the gas, these
forming exoplanets cause deviations from Keplerian motion, which can be
detected through molecular line observations. Our previous work has shown that
machine learning can correctly determine if a planet is present in these
disks. Using our machine learning models, we identify strong, localized non-
Keplerian motion within the disk HD 142666. Subsequent hydrodynamics
simulations of a system with a 5 Jupiter-mass planet at 75 au recreates the
kinematic structure. By currently established standards in the field, we
conclude that HD 142666 hosts a planet. This work represents a first step
towards using machine learning to identify previously overlooked non-Keplerian
features in protoplanetary disks.

Hydrodynamics — Radiative transfer — Accretion disks — Methods: numerical —
Catalogs — Planets and satellites: formation

##  1 Introduction

Protoplanetary accretion disks are the sites of planet formation. The newest
generation of telescopes, such as the Atacama Large Millimeter/submillimeter
Array (ALMA), have unprecedented capabilities for observing protoplanetary
disks. For the first time, we can not only resolve disks themselves, but also
quantify the motion of the dust and gas within them. Disks display a striking
variety of structures such as rings (ALMA Partnership et al., 2015; Dipierro
et al., 2018), likely caused by dust trapping due to forming planets (Pinilla
et al., 2012; Dipierro et al., 2015), and spirals (Pérez et al., 2016), which
may be caused by forming planets (e.g. Dong et al., 2015b) or another
mechanism such as gravitational instability (e.g. Dong et al., 2015a; Hall et
al., 2018; Meru et al., 2017). This new information has greatly advanced our
understanding of the processes underlying the formation and evolution of
planetary systems.

Planets and physical processes, such as gravitational instability, influence
the motion within the disk. This causes the material to deviate from simple
Keplerian motion. Comparing the observed motion against purely Keplerian
motion provides information on the bodies and processes present in the disk
(Hall et al., 2020; Paneque-Carreño et al., 2021; Longarini et al., 2021;
Pinte et al., 2022; Bae et al., 2022; Terry et al., 2022b). Non-Keplerian
motion has been used to uncover a variety of structures, including localized
perturbations associated with gaps and planets (Teague et al., 2018; Pinte et
al., 2018, 2019, 2020) as predicted by Perez et al. (2015).

Kinematic analysis is limited by our ability to accurately identify non-
Keplerian motion. The deviations can be small and frequently occur in noisy
images. It is therefore not only difficult and slow to identify them, but
there is also the strong possibility of overlooking their occurrence. Any
signature that is overlooked is a missed opportunity to detect either a
forming planet or some other process, such as the GI-Wiggle indicative of
gravitational instability (Hall et al., 2020) or the vertical shear
instability (Barraza-Alfaro et al., 2021).

Machine learning (ML) provides a useful tool for this task. ML has quickly
become ubiquitous in both society and the sciences, everything from self-
driving cars (Bojarski et al., 2016) to medicine (Parmar et al., 2015). Recent
efforts in astronomy have made it clear that machine learning is a powerful
method even with simulated training data (Jo & Kim, 2019; Möller & de
Boissière, 2020; Alexander et al., 2020). Machine learning, and in particular
computer vision, excels at the analysis of images (Voulodimos et al., 2018).
In some cases, it has even been shown to outperform humans (Zhou et al.,
2021). It is therefore naturally suited for application to the noisy datasets
in observational astronomy.

Using ML models developed in a previous work (Terry et al., 2022a), we
identify a strong and localized deviation from Keplerian motion in HD 142666.
Using the current widely accepted field standard method (Teague et al., 2018;
Pinte et al., 2018, 2019), we perform smoothed particle hydrodynamic (SPH)
simulations to recreate the kinematic structure of the disk. The agreement is
significant when a 5 MJ planet is included at 75 au. We conclude that HD
142666 hosts a planet.

The paper is arranged as follows: Section 2 describes the models and
simulations used. Section 3 shows the results of applying the models and
simulating the system. Section 4 gives our conclusions.

##  2 Methods

###  2.1 Machine Learning

We use the ML models described in Terry et al. (2022a) and describe them here
for completeness. We use two different architectures: EfficientNetV2 (Tan &
Le, 2021) and RegNet (Xu et al., 2022). All models were made using PyTorch
(Paszke et al., 2019), albeit with significant modifications to the default
models and hyperparameters. We denote these models as EN47, EN61, EN75, RN47,
RN61, and RN75. Table 1 gives performance metrics for the models: model
accuracy at 50% and 95% decision thresholds and the area under the receiver
operating characteristic curve (AUC).

The models were trained using synthetic observations from the MCFOST (Pinte et
al., 2006, 2009) radiative transfer code. MCFOST inputs were drawn from 1000
PHANTOM (Price et al., 2018) SPH simulations of systems with and without
planets (Terry et al., 2022a). Each MCFOST calculation outputs a position-
position-velocity cube from 13CO transition lines ($J=2\rightarrow 1$ and
$J=3\rightarrow 2$). The cubes were convolved spatially and spectrally and
noise was added in order to replicate current observational capabilities.

The model inputs (i.e. radiative transfer outputs) are images of dimension
$C\times H\times W$, where $C$ is the number of input channels, $H$ is the
height of the image, and $W$ is the width of the image (here, $H=W=600$
pixels). A typical grayscale or RGB image will have $C$=1 or 3, respectively.
We instead input an entire position-position-velocity cube.

Observations vary significantly in the number of channels that cover the disk,
but the typical range is between $\approx$40-100. To address this, we train
three different implementations of each model, which gives us a total of six
models. The difference between each implementation is the number of input
velocity channels ($C$=47, 61, or 75).

Each model outputs a two-component vector such that the sum of the components
is 1, i.e. it has undergone softmax activation (Goodfellow et al., 2016). This
can be interpreted as the probability that the given input belongs to a
certain class, i.e. planet- vs no-planet class. The models also output images
of their internal activation structure, which we consider to be the more
important output in this context. While the models were not trained to
pinpoint the locations of planets —a job more suited for semantic image
segmentation (Minaee et al., 2022) —the activation structure can inform us
which regions the model finds important when making its classification
decision. Terry et al. (2022a) found that the activations were able to
highlight velocity channels with non-Keplerian motion in systems that host
planet(s).

To this end, we apply our previously trained models to ALMA data of the HD
142666 system. We inspect the softmax values and activation structures to gain
insight into whether a planet might be in the system and, if so, where its
signature is the strongest.

###  2.2 Observational Data

The HD 142666 data was taken from the DSHARP catalogue (Andrews et al., 2018;
Huang et al., 2018). Data includes 12CO line emission ($J=2\rightarrow 1$) and
1.25 mm continuum images. The system was imaged with a beam with FWHM of
$77\times 61$ mas ($\approx$$11\times 9$ au) with an RMS noise of 1.3 mJy
beam-1; channels have a 0.35 km s-1 resolution (Andrews et al., 2018). Figure
1 shows selected velocity channels overlaid on the continuum.

The image was cropped to focus on the disk, and a subset of velocity channels
was used. The channels were reshaped to $600\times 600$ pixels and normalized
such that all pixel values were between 0 and 1.

###  2.3 Hydrodynamical Simulations

We run a suite of SPH simulations using PHANTOM, varying the mass of the
embedded planet between 1 and 5 MJ. For each simulation, we create channel
maps using MCFOST in the same way that the original training data was made.
The kink is approximately 75 au from the center of the disk, so we place a
planet at this distance.

System parameters are taken from Rubinstein et al. (2018); Andrews et al.
(2018); Huang et al. (2018). The stellar mass, temperature, and radius are 2.0
M⊙, 7500 K, and 2.2 R⊙, respectively. The disk has a mass of 0.0533 M⊙, an
inner radius of 1.3 au, and an outer radius of 150 au. The system is inclined
at 62 degrees with a position angle of 162 and an azimuth of 72 degrees. It is
located 148 pc from Earth.

The SPH outputs are used to create line emission maps to mimic ALMA
capabilities. These calculations are done using the MCFOST radiative transfer
code (Pinte et al., 2006, 2009). Each calculation uses $10^{8}$ photon packets
and includes carbon/silicate dust (Draine & Lee, 1984) with a dust-to-gas
ratio of 1:100. The resulting outputs were convolved spatially and spectrally
to match the observed line emission resolution.

Value | EN47 | EN61 | EN75 | RN47 | RN61 | RN75
---|---|---|---|---|---|---
Accuracy at 50% cutoff (%) | $97\pm 0.5$ | $97\pm 0.5$ | $93\pm 0.7$ | $78\pm 1.1$ | $98\pm 0.4$ | $95\pm 0.6$
Accuracy at 95% cutoff (%) | $96\pm 0.5$ | $94\pm 0.5$ | $88\pm 0.9$ | $65\pm 1.3$ | $96\pm 0.6$ | $92\pm 0.7$
AUC | $0.99\pm 0.002$ | $0.99\pm 0.003$ | $0.98\pm 0.003$ | $0.86\pm 0.010$ | $>0.99\pm 0.001$ | $0.98\pm 0.032$
Table 1: Model performance metrics from Terry et al. (2022a).

##  3 Results and Discussion

Figure 1: Line emission overlaid on continuum. Left: $\Delta v=-1.4$ km/s
channel. Middle: $\Delta v=-1.75$ channel. Right: $\Delta v=-2.1$ channel. The
continuum beam is in magenta, and the line emission beam is in cyan. Figure 2:
HD 142666 structure (12CO: $J=2\rightarrow 1$) and activations. Upper left:
$\Delta v=-1.4$ km/s channel with kink circled in white. Upper middle: $\Delta
v=-1.75$ channel with kink circled in white. Upper right: $\Delta v=-2.1$
channel. Bottom row: selected mean-subtracted activations that roughly
correspond to the channels in the upper row. Activations are from three
different models (EN61, EN47, and RN61, respectively). Line emission beams are
the cyan ellipses in the lower right of the upper row panels. Figure 3: HD
142666 simulation results. Upper left: $\Delta v=-2.3$ km/s channel from the
simulation (convolved beam in lower right). Upper right: $\Delta v=-2.0$ km/s
channel from the simulation (convolved beam in lower right). Lower left:
observed continuum overlaid with contours of simulated $\Delta v=-2.3$ km/s
channel. Lower right: observed continuum overlaid with simulated (cyan) and
observed (white) channels. Continuum beam is in magenta, and the line emission
(simulated and observed) beam is in cyan. The system includes a 5 MJ planet at
75 au. The simulated channels have the continuum and background subtracted for
clarity. The planet’s location is indicated with an x.

Figure 2 shows that HD 142666 has a strong, localized kink that is detected by
the ML models. The kink is particularly visible in the upper middle ($\Delta
v=-1.75$ km/s) channel. The lower row shows activation structures that roughly
correspond to the above channels.

The average softmax value is over 0.84, which means that the models predict
the probability that the input for HD 142666 contains a planet to be over 84%.
This prompts further scrutiny of the activations, which we use to determine
the most probable channel that contains the kink.

The strength and localization of the newly identified kink are reminiscent of
the kinks in HD 163296 and HD 97048. As with HD 163296, the kink in the gas is
outside of the radial extent of the continuum disk. Both of these disks were
found to host planets after SPH simulations containing a planet recreated the
kinematic structure observed in CO observations (Teague et al., 2018; Pinte et
al., 2018, 2019). We apply this same method to HD 142666 to demonstrate that
the kink identified by our models is consistent with kinks identified by
conventional means in HD 163296 and HD 97048.

We found that a simulation of a protoplanetary disk with a 5 MJ planet most
accurately reproduced the observation. Figure 3 shows the results. A localized
kink in the vicinity of the planet is clear in the upper left panel Figure 3
($\Delta v=-2.3$ km/s). This kink is visible to a lesser extent in the $\Delta
v=-2.0$ km/s in the upper right panel of Figure 3, which is also the case in
Figure 1. There is strong agreement between this feature and the non-Keplerian
channel identified by our models: both display a kink of approximately the
same shape and size at approximately the same radial location. This can be
seen in the lower left and right panels of Figure 3. Note that the simulation
and observation do not display the strongest kink in the same velocity
channel. This is simply a relic of the finite temporal resolution of the
simulation, which makes it extremely unlikely that the simulation will be
saved when the planet is exactly coincident with the observation. The temporal
resolution of the simulation was increased to mitigate this effect, but it
persists to some extent.

We conclude that HD 142666 hosts a planet.

We note that our conclusion is confirmed using the same methods described by
Teague et al. (2018); Pinte et al. (2018, 2019). However, what is new about
our approach is that the non-Keplerian motion was first identified by ML
models, highlighting a protoplanet candidate that had previously been missed
upon visual analysis. Verification of the evidence is still done using the
same methodology as previous works (Pinte et al., 2018, 2019). We strongly
advocate that this should always be done for any potential discovery.

###  3.1 Future Work and Limitations

This work shows that machine learning can effectively identify non-Keplerian
motion even if it is missed by humans. However, our work can be improved upon.
The primary limitation is the fact that localising non-Keplerian motion was
not the explicit goal of these models when they were trained. Their purpose
was classification without any attempt of segmentation or object detection.
Models specifically designed to pinpoint deviations would likely be more
effective. Rather than inspecting activation structures —of which there can be
hundreds —the model would directly output a prediction of the location. This
would be a more precise and straightforward method to detect the non-Keplerian
signature, but it would not remove the need to perform follow-up simulations.
We intend to explore this possibility in future works.

Networks such as PGNets (Zhang et al., 2022) and DPNNet-2.0 (Auddy et al.,
2021) offer a potentially fruitful route that would increase the accuracy and
speed of the analysis of disks and channels highlighted by our models. These
networks are designed to infer planetary mass from continuum images. One could
use our models to determine if it is likely that a disk hosts a planet and, if
so, feed the corresponding continuum images into the secondary networks. The
predicted planet mass could then be used as a starting point for followup
simulations rather than simply starting an uninformed parameter sweep. This
would speed up the verification step of the procedure. Such a pipeline could
be useful to explore in future works.

##  4 Conclusion

We have applied ML models created by Terry et al. (2022a) to the DSHARP data
of HD 142666. All models strongly predict the presence of at least one planet.
The activation structures highlight a strong, unreported, and localized kink.
An SPH simulation with a 5 MJ planet at 75 au is able to recreate the newly
identified kinematic structure. By the previously established benchmarks and
methods for kinematic planet detection, we conclude that HD 142666 hosts a
planet.

This work demonstrates the utility of applying machine learning to the
analysis of protoplanetary disks. By highlighting non-Keplerian features in
the disk, ML models are able to guide planet-detection efforts. The signatures
of the planet were previously overlooked by human analysts, and the
traditional analysis was only performed because of the information given by
the models. We anticipate that this method can identify new non-Keplerian
features in both existing and future protoplanetary observations.

##  5 Acknowledgements

This paper makes use of the following ALMA data: ADS/JAO.ALMA #2016.1.00484.L.
ALMA is a partnership of ESO (representing its member states), NSF (USA) and
NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ. J.T. was a participant in
the 2022 Machine Learning for Science (ML4SCI) Google Summer of Code program.
S.G. was supported in part by the National Science Foundation Award No.
2108645. This study was supported in part by resources and technical expertise
from the Georgia Advanced Computing Resource Center, a partnership between the
University of Georgia’s Office of the Vice President for Research and Office
of the Vice President for Information Technology.

## References

  * Alexander et al. (2020) Alexander, S., Gleyzer, S., McDonough, E., Toomey, M. W., & Usai, E. 2020, ApJ, 893, 15, doi: 10.3847/1538-4357/ab7925
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3, doi: 10.1088/2041-8205/808/1/L3
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Auddy et al. (2021) Auddy, S., Dey, R., Lin, M.-K., & Hall, C. 2021, ApJ, 920, 3, doi: 10.3847/1538-4357/ac1518
  * Bae et al. (2022) Bae, J., Teague, R., Andrews, S. M., et al. 2022, ApJ, 934, L20, doi: 10.3847/2041-8213/ac7fa3
  * Barraza-Alfaro et al. (2021) Barraza-Alfaro, M., Flock, M., Marino, S., & Pérez, S. 2021, A&A, 653, A113, doi: 10.1051/0004-6361/202140535
  * Bojarski et al. (2016) Bojarski, M., Del Testa, D., Dworakowski, D., et al. 2016, arXiv e-prints, arXiv:1604.07316.  https://arxiv.org/abs/1604.07316
  * Dipierro et al. (2015) Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73, doi: 10.1093/mnrasl/slv105
  * Dipierro et al. (2018) Dipierro, G., Ricci, L., Pérez, L., et al. 2018, MNRAS, 475, 5296, doi: 10.1093/mnras/sty181
  * Dong et al. (2015a) Dong, R., Hall, C., Rice, K., & Chiang, E. 2015a, ApJ, 812, L32, doi: 10.1088/2041-8205/812/2/L32
  * Dong et al. (2015b) Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015b, ApJ, 809, L5, doi: 10.1088/2041-8205/809/1/L5
  * Draine & Lee (1984) Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89, doi: 10.1086/162480
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., & Courville, A. 2016, Deep Learning (MIT Press)
  * Hall et al. (2018) Hall, C., Rice, K., Dipierro, G., et al. 2018, MNRAS, 477, 1004, doi: 10.1093/mnras/sty550
  * Hall et al. (2020) Hall, C., Dong, R., Teague, R., et al. 2020, ApJ, 904, 148, doi: 10.3847/1538-4357/abac17
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42, doi: 10.3847/2041-8213/aaf740
  * Jo & Kim (2019) Jo, Y., & Kim, J.-h. 2019, MNRAS, 489, 3565, doi: 10.1093/mnras/stz2304
  * Longarini et al. (2021) Longarini, C., Lodato, G., Toci, C., et al. 2021, ApJ, 920, L41, doi: 10.3847/2041-8213/ac2df6
  * Meru et al. (2017) Meru, F., Juhász, A., Ilee, J. D., et al. 2017, ApJ, 839, L24, doi: 10.3847/2041-8213/aa6837
  * Minaee et al. (2022) Minaee, S., Boykov, Y., Porikli, F., et al. 2022, IEEE Transactions on Pattern Analysis and Machine Intelligence, 44, 3523, doi: 10.1109/TPAMI.2021.3059968
  * Möller & de Boissière (2020) Möller, A., & de Boissière, T. 2020, MNRAS, 491, 4277, doi: 10.1093/mnras/stz3312
  * Paneque-Carreño et al. (2021) Paneque-Carreño, T., Pérez, L. M., Benisty, M., et al. 2021, ApJ, 914, 88, doi: 10.3847/1538-4357/abf243
  * Parmar et al. (2015) Parmar, C., Grossmann, P., Bussink, J., Lambin, P., & Aerts, H. J. W. L. 2015, Scientific Reports, 5, 13087, doi: 10.1038/srep13087
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., et al. 2019, in Advances in Neural Information Processing Systems, ed. H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, & R. Garnett, Vol. 32 (Curran Associates, Inc.), doi: 10.48550/arXiv.1912.01703
  * Pérez et al. (2016) Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519, doi: 10.1126/science.aaf8296
  * Perez et al. (2015) Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJ, 811, L5, doi: 10.1088/2041-8205/811/1/L5
  * Pinilla et al. (2012) Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114, doi: 10.1051/0004-6361/201118204
  * Pinte et al. (2009) Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967, doi: 10.1051/0004-6361/200811555
  * Pinte et al. (2006) Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797, doi: 10.1051/0004-6361:20053275
  * Pinte et al. (2022) Pinte, C., Teague, R., Flaherty, K., et al. 2022, arXiv e-prints, arXiv:2203.09528.  https://arxiv.org/abs/2203.09528
  * Pinte et al. (2018) Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13, doi: 10.3847/2041-8213/aac6dc
  * Pinte et al. (2019) Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nature Astronomy, 3, 1109, doi: 10.1038/s41550-019-0852-6
  * Pinte et al. (2020) Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9, doi: 10.3847/2041-8213/ab6dda
  * Price et al. (2018) Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, PASA, 35, e031, doi: 10.1017/pasa.2018.25
  * Rubinstein et al. (2018) Rubinstein, A. E., Macías, E., Espaillat, C. C., et al. 2018, ApJ, 860, 7, doi: 10.3847/1538-4357/aabfba
  * Tan & Le (2021) Tan, M., & Le, Q. 2021, in Proceedings of Machine Learning Research, Vol. 139, Proceedings of the 38th International Conference on Machine Learning, ed. M. Meila & T. Zhang (PMLR), 10096–10106, doi: 10.48550/arXiv.2104.00298
  * Teague et al. (2018) Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12, doi: 10.3847/2041-8213/aac6d7
  * Terry et al. (2022a) Terry, J. P., Hall, C., Abreau, S., & Gleyzer, S. 2022a, ApJ, 941, 192, doi: 10.3847/1538-4357/aca477
  * Terry et al. (2022b) Terry, J. P., Hall, C., Longarini, C., et al. 2022b, MNRAS, 510, 1671, doi: 10.1093/mnras/stab3513
  * Voulodimos et al. (2018) Voulodimos, A., Doulamis, N., Doulamis, A., Protopapadakis, E., & Andina, D. 2018, Neuroscience, 7068349, doi: 10.1155/2018/7068349
  * Xu et al. (2022) Xu, J., Pan, Y., Pan, X., et al. 2022, IEEE Transactions on Neural Networks and Learning Systems, 1, doi: 10.1109/TNNLS.2022.3158966
  * Zhang et al. (2022) Zhang, S., Zhu, Z., & Kang, M. 2022, MNRAS, 510, 4473, doi: 10.1093/mnras/stab3502
  * Zhou et al. (2021) Zhou, W., Yang, Y., Yu, C., et al. 2021, Nature Communications, 12, 1259, doi: 10.1038/s41467-021-21466-z