# A Fast Template Periodogram for Detecting Non-Sinusoidal Fixed-Shape Signals
in Irregularly Sampled Time Series

J. Hoffman11affiliation: Xaxis, LLC, 3 World Trade Center, 175 Greenwich
Street, 30th Floor, New York, NY 10007 , J. Vanderplas22affiliation: eScience
Institute, University of Washington, Seattle, WA 98195 , J. D.
Hartman33affiliation: Department of Astrophysical Sciences, Princeton
University, Princeton NJ 08540 , G. Á Bakos33affiliation: Department of
Astrophysical Sciences, Princeton University, Princeton NJ 08540
44affiliation: Institute for Advanced Study, 1 Einsten Drive, Princeton, NJ
<EMAIL_ADDRESS>
###### Abstract

Astrophysical time series often contain periodic signals. The large and
growing volume of time series data from photometric surveys demands
computationally efficient methods for detecting and characterizing such
signals. The most efficient algorithms available for this purpose are those
that exploit the $\mathcal{O}(N\log N)$ scaling of the Fast Fourier Transform
(FFT). However, these methods are not optimal for non-sinusoidal signal
shapes. Template fits (or periodic matched filters) optimize sensitivity for
_a priori_ known signal shapes but at a significant computational cost.
Current implementations of template periodograms scale as
$\mathcal{O}(N_{f}N_{\rm obs})$, where $N_{f}$ is the number of trial
frequencies and $N_{\rm obs}$ is the number of lightcurve observations, and
due to non-convexity, they do not guarantee the best fit at each trial
frequency, which can lead to spurious results. In this work, we present a non-
linear extension of the Lomb-Scargle periodogram to obtain a template-fitting
algorithm that is both accurate (globally optimal solutions are obtained
except in pathological cases) and computationally efficient (scaling as
$\mathcal{O}(N_{f}\log N_{f})$ for a given template). The non-linear
optimization of the template fit at each frequency is recast as a polynomial
zero-finding problem, where the coefficients of the polynomial can be computed
efficiently with the non-equispaced fast Fourier transform. We show that our
method, which uses truncated Fourier series to approximate templates, is an
order of magnitude faster than existing algorithms for small problems
($N\lesssim 10$ observations) and 2 orders of magnitude faster for long base-
line time series with $N_{\rm obs}\gtrsim 10^{4}$ observations. An open-source
implementation of the fast template periodogram is available at
github.com/PrincetonUniversity/FastTemplatePeriodogram.

##  1\. Introduction

Astronomical systems exhibit a wide range of time-dependent variability. By
measuring and characterizing this variability, astronomers are able to infer a
variety of important astrophysical properties about the underlying system.

Periodic signals in noisy astronomical timeseries can be detected with a
number of techniques, including Gaussian process regression (Foreman-Mackey et
al., 2017; Rasmussen & Williams, 2005), least-squares spectral analysis (Lomb,
1976; Scargle, 1982; Barning, 1963; Vaníček, 1971), and information-theoretic
methods (Graham et al., 2013a; Huijse et al., 2012; Cincotta et al., 1995).
For an empirical comparison of some of these techniques applied to several
astronomical survey datasets, see Graham et al. (2013b).

For stationary periodic signals, least-squares spectral analysis — also known
as the Lomb-Scargle (LS) periodogram (Lomb, 1976; Scargle, 1982; Barning,
1963; Vaníček, 1971) — is perhaps the most sensitive and computationally
efficient method of detection. The LS periodogram can be made to scale as
$\mathcal{O}(N_{f}\log N_{f})$, where $N_{f}$ is the number of trial
frequencies, by utilizing the non-equispaced fast Fourier transform (Keiner et
al., 2009; Dutt & Rokhlin, 1993, NFFT) to evaluate frequency-dependent sums,
or by “extirpolating” irregularly spaced observations to a regular grid with
Lagrange polynomials (Press & Rybicki, 1989).

The LS periodogram fits the following model to a set of observations:

$\hat{y}_{\rm LS}(t|\theta,\omega)=\theta_{0}\cos{\omega
t}+\theta_{1}\sin{\omega t}.$ (1)

where $\omega=2\pi/P$ is the (angular) frequency of the underlying signal, and
$\theta_{0}$ and $\theta_{1}$ are the amplitudes of the signal. When the data
$y_{i}$ is composed of a sinusoidal component and a white noise component
(i.e., when the measurement uncertainties are uncorrelated and Gaussian), the
LS periodogram provides a maximum likelihood estimate for the model parameters
($\omega,\theta_{0},$ and $\theta_{1}$).

The LS “power” $P_{LS}(\omega)$ has several definitions in the literature
(Zechmeister & Kürster, 2009), but we adopt the following definition
throughout the paper:

$P(\omega)=\frac{\chi^{2}_{0}-\chi^{2}(\omega)}{\chi^{2}_{0}}$ (2)

Where $\chi^{2}_{0}$ is the weighted sum of squared residuals for a constant
fit:

$\chi^{2}_{0}=\sum_{i}w_{i}(y_{i}-\bar{y})^{2}$ (3)

where $\bar{y}=\sum_{i}w_{i}y_{i}$ is the weighted mean of the observations,
and $w_{i}\propto\sigma_{i}^{-2}$ are the normalized weights for each
observation ($\sum_{i}w_{i}=1$), and $\chi^{2}(\omega)$ is the weighted sum of
squared residuals for the best-fit model $\hat{y}$ at a given trial frequency
$\omega=2\pi/P$ where $P$ is the period:

$\chi^{2}(\omega)=\min_{\theta}\sum_{i}w_{i}(y_{i}-\hat{y}(t_{i}|\omega,\theta))^{2}$
(4)

###  1.1. Bayesian interpretation

We note that, while this formalism captures many data and modeling scenarios,
it is not completely general. For example, correlated uncertainties are not
handled here. A more general Bayesian treatment of periodic models in
astronomical timeseries is better handled by expressing a posterior over the
model parameters.

Assuming a Gaussian likelihood for the observations
$y_{i}\sim\mathcal{N}(\hat{y}(t_{i}|\omega,\theta),\sigma^{2}_{i})$, and
uniform priors on both the frequency parameter $\omega$ and the non-frequency
parameters $\theta$, the posterior is

$\displaystyle p(\omega,\theta|X)$
$\displaystyle=\frac{p(X|\omega,\theta)p(\omega,\theta)}{p(X)}$ (5)
$\displaystyle\propto
p(\omega,\theta)\prod_{i=1}^{N_{\mathrm{obs}}}\mathcal{N}(\hat{y}(t_{i}|\omega,\theta),\sigma^{2}_{i})$
(6)

where we use $X=\\{(t_{i},y_{i},\sigma_{i})|0<i<N_{\mathrm{obs}}\\}$ as
shorthand for the lightcurve observations. The logarithm of the posterior is

$\displaystyle\log p(\omega,\theta|X)$
$\displaystyle=-\frac{1}{2}\sum_{i=1}^{N_{\mathrm{obs}}}\left(\frac{y_{i}-\hat{y}(t_{i}|\omega,\theta)}{\sigma_{i}}\right)^{2}+\mathrm{const.}$
(7)
$\displaystyle=-\frac{W}{2}\sum_{i=1}^{N_{\mathrm{obs}}}w_{i}(y_{i}-\hat{y}(t_{i}|\omega,\theta))^{2}+\mathrm{const.}$
(8)

where $W=\sum_{i}\sigma_{i}^{-2}$. Thus,

$\displaystyle\chi^{2}(\omega)$
$\displaystyle=\min_{\theta}\sum_{i}w_{i}(y_{i}-\hat{y}(t_{i}|\omega,\theta))^{2}$
(9) $\displaystyle=\min_{\theta}\left(-\frac{2}{W}\log
p(\omega,\theta|X)+\mathrm{const.}\right)$ (10)
$\displaystyle=-\frac{2}{W}\max_{\theta}\log
p(\omega,\theta|X)+\mathrm{const.}$ (11)

and therefore, since $P(\omega)=1-\chi^{2}(\omega)/\chi_{0}^{2}$,

$P(\omega)=\frac{2}{W\chi_{0}^{2}}\max_{\theta}p(\omega,\theta|X)+\mathrm{const.}$
(12)

The Lomb-Scargle power is a linear transformation of the maximum of the log
posterior over the non-frequency parameters, with the frequency parameter held
fixed. Thus, choosing the frequency that maximizes the periodogram value
corresponds to finding a MAP estimate of the frequency parameter.

A MAP interpretation is more general, and is applicable to scenarios not
considered in this paper (e.g. correlated uncertainties, multi-dimensional
timeseries, etc.), since all of these problems are amenable to MAP estimation
of their model parameters. However, in order to maintain consistency with
notation in the Lomb-Scargle literature (e.g. Zechmeister & Kürster, 2009;
VanderPlas, 2018), we keep our definition of the periodogram restricted as
above and only consider one-dimensional timeseries with heteroscedastic but
uncorrelated uncertainties.

###  1.2. Extending Lomb-Scargle

The LS periodogram has numerous extensions to account for, e.g., biased
estimates of the mean brightness (Zechmeister & Kürster, 2009), non-sinusoidal
signals (Schwarzenberg-Czerny, 1996; Palmer, 2009), multi-band observations
(VanderPlas & Ivezić, 2015), and to mitigate overfitting of more flexible
models via regularization (VanderPlas & Ivezić, 2015). For a detailed review
of the LS periodogram and its extensions, see VanderPlas (2018).

The multi-harmonic LS periodogram (Bretthorst et al., 1988; Schwarzenberg-
Czerny, 1996; Palmer, 2009, MHGLS), provides a more flexible model by adding
harmonic components to the fit:

$\hat{y}_{\rm
MHLS}(t|\theta,\omega)=\theta_{0}+\sum_{n=1}^{H}\theta_{2n}\cos{n\omega
t}+\theta_{2n-1}\sin{n\omega t}.$ (13)

Additional harmonics are important for modeling signals that, while stationary
and periodic, are non-sinusoidal (e.g. RR Lyrae, eclipsing binaries, etc.).

However, the MHLS periodogram contains $2H+1$ free parameters while the
original LS periodogram contains only $2$ ($3$ if the mean brightness is
considered a free parameter). Including higher order harmonics adds model
complexity, which can degrade the sensitivity of the MHLS periodogram to
sinusoidal or approximately sinusoidal signals.

Tikhonov regularization (or $L_{2}$ regularization), is one tool for
mitigating overfitting of the higher order harmonics. However, adding an
$L_{2}$ regularization term to the Fourier amplitudes adds bias to the model,
and that bias should be compared to the value added by the decrease in model
variance (VanderPlas & Ivezić, 2015).

###  1.3. Computational scaling

The LS periodogram naively scales as $\mathcal{O}(N_{f}N_{\rm obs})$ where
$N_{f}$ is the number of trial frequencies and $N_{\rm obs}$ is the number of
observations. However, the limiting computations of the LS periodogram involve
sums of trigonometric functions over the observations. When the observations
are regularly sampled, the fast Fourier transform (FFT) (Cooley & Tukey, 1965)
can evaluate such sums efficiently and the LS periodogram scales as
$\mathcal{O}(N_{f}\log N_{f})$.

When the data is not regularly sampled, as is the case for most astronomical
time series, the LS periodogram can be evaluated quickly in one of two popular
ways. The first, by Press & Rybicki (1989) involves“extirpolating” irregularly
sampled data onto a regularly sampled mesh, and then performing FFTs to
evaluate the necessary sums. The second, as pointed out in Leroy (2012), is to
use the non-equispaced FFT (Keiner et al., 2009; Dutt & Rokhlin, 1993, NFFT)
to evaluate the sums; this provides roughly an order of magnitude speedup over
the Press & Rybicki (1989) algorithm, and both algorithms scale as
$\mathcal{O}(N_{f}\log N_{f})$ (Leroy, 2012).

There is a growing population of alternative methods for detecting periodic
signals in astrophysical data. Some of these methods can reliably outperform
the LS periodogram, especially for non-sinusoidal signal shapes (see Graham et
al. (2013b) for a recent empirical review of period finding algorithms).
However, a key advantage the LS periodogram and its extensions is speed.
Virtually all other “phase-folding” methods scale as $\mathcal{O}(N_{\rm
obs}\times N_{f})$, where $N_{\rm obs}$ is the number of observations and
$N_{f}$ is the number of trial frequencies, while the Lomb-Scargle periodogram
scales as $\mathcal{O}(N_{f}\log N_{f})$. The virtual independence of Lomb-
Scargle’s computation time with respect to the number of observations
(assuming $N_{f}\gtrsim N_{\rm obs}$) is especially valuable for lightcurves
with $N_{\rm obs}\gg\log N_{f}\sim 50$.

Algorithmic efficiency will become increasingly important as the volume of
data produced by astronomical observatories continues to grow larger. The
HATNet survey (Bakos et al., 2004), for example, has already made
$\mathcal{O}(10^{4})$ observations of $\mathcal{O}(10^{6}-10^{7})$ stars. The
Gaia telescope (Gaia Collaboration et al., 2016) is set to produce
$\mathcal{O}(10-100)$ observations of $\mathcal{O}(10^{9})$ stars. The Large
Synoptic Survey Telescope (LSST; LSST Science Collaboration et al. (2009))
will make $\mathcal{O}(10^{2}-10^{3})$ observations of $\mathcal{O}(10^{10})$
stars during its operation starting in 2023.

###  1.4. Template periodograms

When the shape of a stationary periodic signal is known a priori, then the
number of degrees of freedom is the same as the original LS periodogram (with
a floating mean component):

$\hat{y}(t|\theta,\omega)=\theta_{0}+\theta_{1}\mathbf{M}(\omega
t-\theta_{2}),$ (14)

where $\mathbf{M}:[0,2\pi)\rightarrow\mathcal{R}$ is a predefined periodic
template. We refer to the periodogram corresponding to this model as the
“template periodogram.”

As is the case for the LS periodogram, the template periodogram is equivalent
to a maximum-likelihood estimate of the model parameters under the assumption
that measurement uncertainties are Gaussian and uncorrelated (i.e. white
noise).

This paper develops new extensions of least-squares spectral analysis for
arbitrary signal shapes. For non-periodic signals this method is known as
matched filter analysis, and can be extended to search for periodic signals
by, e.g., phase folding the data at different trial periods.

An analysis by Sesar et al. (2016) found that template fitting significantly
improved period and amplitude estimation for RR Lyrae in Pan-STARRS DR1
photometry (Chambers et al., 2016). Since the signal shapes for RR Lyrae in
various bandpasses are known _a priori_ (see Sesar et al. (2010)), template
fitting provides an optimal estimate of amplitude and period, given that the
object is indeed an RR Lyrae star well modeled by at least one of the
templates. Templates were especially crucial for Pan-STARRS data, since there
are typically only 35 observations per source over 5 bands (Hernitschek et
al., 2016), not enough to obtain accurate amplitudes empirically by phase-
folding. By including domain knowledge (i.e. knowledge of what RR Lyrae
lightcurves look like), template fitting allows for accurate inferences of
amplitude even for undersampled lightcurves.

However, the improved accuracy comes at substantial computational cost: the
template fitting procedure took 30 minutes per CPU per object, and Sesar et
al. (2016) were forced to limit the number of fitted lightcurves ($\lesssim
1000$) in order to keep the computational costs to a reasonable level. Several
cuts were made before the template fitting step to reduce the more than 1
million Pan-STARRS DR1 objects to a small enough number, and each of these
steps removes a small portion of RR Lyrae from the sample. Though this number
was reported by Sesar et al. (2016) to be small ($\lesssim 2\%$), it may be
possible to further improve the completeness of the final sample by applying
template fits to a larger number of objects, which would require either more
computational resources, more time, or, ideally, a more efficient template
fitting procedure.

The paper is organized as follows. Section 2 poses the problem of template
fitting in the language of least squares spectral analysis and derives the
fast template periodogram. Section 3 describes a freely available
implementation of the new template periodogram. Section 4 summarizes our
results, addresses caveats, and discusses possible avenues for improving the
efficiency of the current algorithm.

##  2\. Derivations

We define a template $\mathbf{M}$

$\mathbf{M}:[0,2\pi)\rightarrow\mathbb{R},$ (15)

as a mapping between the unit interval and the set of real numbers. We
restrict our discussion to sufficiently smooth templates such that
$\mathbf{M}$ can be adequately described by a truncated Fourier series

$\hat{\mathbf{M}}(\omega t|H)=\sum_{n=1}^{H}\left(c_{n}\cos{n\omega
t}+s_{n}\sin{n\omega t}\right)$ (16)

for some finite $H>0$.

That the $c_{n}$ and $s_{n}$ values are _fixed_ (i.e., they define the
template) is the crucial difference between the template periodogram and the
multi-harmonic Lomb-Scargle (Palmer, 2009; Bretthorst et al., 1988), where
$c_{n}$ and $s_{n}$ are _free parameters_.

We now construct a periodogram for this template. The periodogram assumes that
an observed time series $S=\\{(t_{i},y_{i},\sigma_{i})\\}_{i=1}^{N}$ can be
modeled by a scaled, transposed template that repeats with period
$2\pi/\omega$, i.e.

$y_{i}\approx\hat{y}(\omega
t_{i}|\theta,\mathbf{M})=\theta_{1}\mathbf{M}(\omega
t_{i}-\theta_{2})+\theta_{3},$ (17)

where $\theta=(\theta_{1},\theta_{2},\theta_{3})\in\mathbb{R}^{3}$ is a set of
model parameters.

The optimal parameters are the location of a local minimum of the (weighted)
sum of squared residuals,

$\chi^{2}(\theta,S)\equiv\sum_{i}w_{i}(y_{i}-\hat{y}(\omega
t_{i}|\theta))^{2},$ (18)

and thus the following condition must hold for all three model parameters at
the optimal solution $\theta=\theta_{\rm opt}$:

$\left.\frac{\partial\chi^{2}}{\partial\theta_{j}}\right|_{\theta=\theta_{\rm
opt}}=0~{}~{}\forall\theta_{j}\in\theta.$ (19)

Note that we have implicitly assumed $\chi^{2}(\theta,S)$ is a $C^{1}$
differentiable function of $\theta$, which requires that $\mathbf{M}$ is a
$C^{1}$ differentiable function. Though this assumption could be violated if
we considered a more complete set of templates, (e.g. a box function), our
restriction to truncated Fourier series ensures $C^{1}$ differentiability.

Note that we also implicitly assume $\sigma_{i}>0$ for all $i$ and we will
later assume that the variance of the observations $y_{i}$ is non-zero. If
there are no measurement errors, i.e. $\sigma_{i}=0$ for all $i$, then uniform
weights (setting $\sigma_{i}=1$) should be used. If the variance of the
observations $y$ is zero, the periodogram (as defined in Equation 2) is
undefined for all frequencies. We do not consider the case where
$\sigma_{i}=0$ for some observations $i$ and $\sigma_{j}>0$ for some
observations $j$.

We can derive a system of equations for $\theta_{\rm opt}$ from the condition
given in Equation 19. The explicit condition that must be met for each
parameter $\theta_{j}$ is simplified below, using

$\hat{y}_{i}=\hat{y}(\omega t_{i}|\theta)$ (20)

and

$\partial_{j}\hat{y}_{i}=\left.\frac{\partial\hat{y}(\omega
t|\theta)}{\partial\theta_{j}}\right|_{t=t_{i}}$ (21)

for brevity:

$\begin{split}0&=\left.\frac{\partial\chi^{2}}{\partial\theta_{j}}\right|_{\theta=\theta_{\rm
opt}}\\\
&=-2\sum_{i}w_{i}\left(y_{i}-\hat{y}_{i}\right)(\partial_{j}\hat{y})_{i}\\\
\sum_{i}w_{i}y_{i}(\partial_{j}\hat{y})_{i}&=\sum_{i}w_{i}\hat{y}_{i}(\partial_{j}\hat{y})_{i}.\end{split}$
(22)

The above is a general result that extends to all least squares periodograms.
To simplify derivations, we adopt the following notation:

$\displaystyle\left<X\right>$ $\displaystyle\equiv$
$\displaystyle\sum_{i}w_{i}X_{i}$ (23) $\displaystyle\left<XY\right>$
$\displaystyle\equiv$ $\displaystyle\sum_{i}w_{i}X_{i}Y_{i}$ (24)
$\displaystyle{\rm Cov}(X,Y)$ $\displaystyle\equiv$
$\displaystyle\left<XY\right>-\left<X\right>\left<Y\right>$ (25)
$\displaystyle{\rm Var}(X)$ $\displaystyle\equiv$ $\displaystyle{\rm
Cov}(X,X)$ (26)

In addition, the transposed template
$\mathbf{M}_{\theta_{2}}=\mathbf{M}(\omega t_{i}-\theta_{2})$ can be expressed
as

$\displaystyle\mathbf{M}_{\theta_{2}}(\omega t)$
$\displaystyle=\sum_{n}c_{n}\cos n\left(\omega t-\theta_{2}\right)$ (27)
$\displaystyle\qquad+s_{n}\sin{n\left(\omega t-\theta_{2}\right)}$ (28)
$\displaystyle=\sum_{n}\left(c_{n}\cos{n\theta_{2}}-s_{n}\sin{n\theta_{2}}\right)\cos{n\omega
t}+$ (29)
$\displaystyle\qquad\left(s_{n}\cos{n\theta_{2}}+c_{n}\sin{n\theta_{2}}\right)\sin{n\omega
t}$ (30)
$\displaystyle=\sum_{n}\left(\alpha_{n}e^{in\theta_{2}}+\alpha_{n}^{*}e^{-in\theta_{2}}\right)\cos{n\omega
t}+$ (31)
$\displaystyle\qquad\left(-i\left[\alpha_{n}e^{in\theta_{2}}-\alpha_{n}^{*}e^{-in\theta_{2}}\right]\right)\sin{n\omega
t}$ (32)
$\displaystyle=\sum_{n}\left(\alpha_{n}\psi^{n}+\alpha_{n}^{*}\psi^{-n}\right)\cos{n\omega
t}+$ (33)
$\displaystyle\qquad\left(-i\left[\alpha_{n}\psi^{n}-\alpha_{n}^{*}\psi^{-n}\right]\right)\sin{n\omega
t}$ (34) $\displaystyle=\sum_{n}A_{n}(\psi)\cos{n\omega
t}+B_{n}(\psi)\sin{n\omega t}$ (35)

where $\alpha_{n}=(c_{n}+is_{n})/2$, and $\psi\equiv e^{i\theta_{2}}$ is a
convenient change of variable.

We also define the following terms:

$\displaystyle\widehat{YM}$
$\displaystyle=\left<y\mathbf{M}_{\theta_{2}}\right>$ (36)
$\displaystyle\widehat{MM}$
$\displaystyle=\left<\mathbf{M}_{\theta_{2}}^{2}\right>$ (37)
$\displaystyle\overline{M}$
$\displaystyle=\left<\mathbf{M}_{\theta_{2}}\right>$ (38) $\displaystyle MM$
$\displaystyle={\rm
Var}(\mathbf{M}_{\theta_{2}})=\widehat{MM}-\overline{M}^{2}$ (39)
$\displaystyle YM$ $\displaystyle={\rm
Cov}(\mathbf{M}_{\theta_{2}},y)=\widehat{YM}-\bar{y}\overline{M}$ (40)

For a given phase shift $\theta_{2}$, the optimal amplitude and offset are
obtained from requiring the partial derivatives of the sum of squared
residuals, $\chi^{2}$, to be zero.

Namely, we obtain that

$\displaystyle 0=\frac{\partial\chi^{2}}{\partial\theta_{1}}$
$\displaystyle=2\sum_{i}w_{i}(y_{i}-\hat{y}_{i})\left(-\frac{\partial\hat{y}}{\partial\theta_{1}}\right)_{i}$
(41)
$\displaystyle=\sum_{i}w_{i}(y_{i}-\theta_{1}\mathbf{M}_{\theta_{2}}-\theta_{3})\mathbf{M}_{\theta_{2}}$
(42)
$\displaystyle=\widehat{YM}-\theta_{1}\widehat{MM}-\theta_{3}\overline{M}$
(43)

and

$\displaystyle 0=\frac{\partial\chi^{2}}{\partial\theta_{3}}$
$\displaystyle=2\sum_{i}w_{i}(y_{i}-\hat{y}_{i})\left(-\frac{\partial\hat{y}}{\partial\theta_{3}}\right)_{i}$
(44)
$\displaystyle=\sum_{i}w_{i}(y_{i}-\theta_{1}\mathbf{M}_{\theta_{2}}-\theta_{3})$
(45) $\displaystyle=\bar{y}-\theta_{1}\overline{M}-\theta_{3}$ (46)

This system of equations can then be rewritten as

$\begin{pmatrix}\widehat{MM}&\overline{M}\\\
\overline{M}&1\end{pmatrix}\begin{pmatrix}\theta_{1}\\\
\theta_{3}\end{pmatrix}=\begin{pmatrix}\widehat{YM}\\\ \bar{y}\end{pmatrix}$
(47)

which reduces to

$\displaystyle\begin{pmatrix}\theta_{1}\\\ \theta_{3}\end{pmatrix}$
$\displaystyle=\frac{1}{\widehat{MM}-\overline{M}^{2}}\begin{pmatrix}1&-\overline{M}\\\
-\overline{M}&\widehat{MM}\end{pmatrix}\begin{pmatrix}\widehat{YM}\\\
\bar{y}\end{pmatrix}$ (48)
$\displaystyle=\frac{1}{\widehat{MM}-\overline{M}^{2}}\begin{pmatrix}\widehat{YM}-\bar{y}\overline{M}\\\
\widehat{MM}\bar{y}-\widehat{YM}\overline{M}\end{pmatrix}$ (49)

Letting $MM=\widehat{MM}-\overline{M}^{2}$ and
$YM=\widehat{YM}-\bar{y}\overline{M}$, we have

$\begin{pmatrix}\theta_{1}\\\ \theta_{3}\end{pmatrix}=\begin{pmatrix}YM/MM\\\
\bar{y}-\overline{M}(YM/MM)\end{pmatrix}$ (50)

This means we can rewrite the model
$\hat{y}=\theta_{1}\mathbf{M}_{\theta_{2}}+\theta_{3}$ as

$\hat{y}_{i}=\bar{y}+\left(\frac{YM}{MM}\right)(M_{i}-\overline{M})$ (51)

To obtain an expression for the periodogram, $P=1-\chi^{2}/\chi^{2}_{0}$, we
first compute $\chi^{2}$

$\displaystyle\chi^{2}$ $\displaystyle=\sum_{i}w_{i}(y_{i}-\hat{y}_{i})^{2}$
(52)
$\displaystyle=\sum_{i}w_{i}(y_{i}^{2}-2y_{i}\hat{y}_{i}+\hat{y}_{i}^{2})$
(53) $\displaystyle=YY-2\frac{(YM)^{2}}{MM}+\frac{(YM)^{2}}{MM}$ (54)
$\displaystyle=YY-\frac{(YM)^{2}}{MM}$ (55)

Since, $\chi^{2}_{0}=YY$, we have

$P(\omega)=\frac{(YM)^{2}}{YY\cdot MM}$ (56)

We wish to maximize $P(\omega)$ with respect to the phase shift parameter
$\theta_{2}$,

$\displaystyle\partial_{\theta_{2}}P=0$ $\displaystyle=\frac{YM}{YY\cdot
MM}\left(2\partial_{\theta_{2}}(YM)-\frac{YM}{MM}\partial_{\theta_{2}}(MM)\right)$
(57) $\displaystyle=2MM\partial_{\theta_{2}}(YM)-YM\partial_{\theta_{2}}(MM).$
(58)

The final expression is the non-linear condition that must be satisfied by the
optimal phase shift parameter $\theta_{2}$. However, satisfying Equation 57 is
not _sufficient_ to guarantee that $\theta_{2}$ is optimal. The value of the
periodogram at each $\theta_{2}$ satisfying Equation 57 must be computed, and
the globally optimal solution chosen from this set.

We seek a more explicit form for Equation 57. We derive expressions for $MM$
and $YM$, defining

$\displaystyle CC_{nm}$ $\displaystyle\equiv$ $\displaystyle{\rm
Cov}(\cos{n\omega t},\cos{m\omega t})$ (59) $\displaystyle CS_{nm}$
$\displaystyle\equiv$ $\displaystyle{\rm Cov}(\cos{n\omega t},\sin{m\omega
t})$ (60) $\displaystyle SS_{nm}$ $\displaystyle\equiv$ $\displaystyle{\rm
Cov}(\sin{n\omega t},\sin{m\omega t})$ (61) $\displaystyle YC_{n}$
$\displaystyle\equiv$ $\displaystyle\left<(y-\bar{y})\cos{n\omega t}\right>$
(62) $\displaystyle YS_{n}$ $\displaystyle\equiv$
$\displaystyle\left<(y-\bar{y})\sin{n\omega t}\right>,$ (63)

all of which can be evaluated efficiently using the NFFT.

The autocovariance of the template values $MM$, is given by

$\displaystyle MM$
$\displaystyle\equiv\sum_{i}w_{i}M_{i}^{2}-\left(\sum_{i}w_{i}M_{i}\right)^{2}$
(64) $\displaystyle=\sum_{i}w_{i}\left(\sum_{n}A_{n}\cos{\omega
nt_{i}}+B_{n}\sin{\omega nt_{i}}\right)^{2}$ (65)
$\displaystyle\qquad-\left(\sum_{n}A_{n}C_{n}+B_{n}S_{n}\right)^{2}$ (66)
$\displaystyle=\sum_{n,m}A_{n}A_{m}CC_{nm}$
$\displaystyle\qquad+(A_{n}B_{m}CS_{nm}+B_{n}A_{m}(CS^{T})_{nm})$
$\displaystyle\qquad+B_{n}B_{m}SS_{nm}$ (67)
$\displaystyle=\sum_{n,m}A_{n}A_{m}CC_{nm}+2A_{n}B_{m}CS_{nm}$
$\displaystyle\qquad+B_{n}B_{m}SS_{nm},$ (68)

using

$\sum_{n,m}A_{n}B_{m}CS_{nm}=\sum_{n,m}A_{m}B_{n}CS_{mn}.$ (69)

We also derive the products $A_{n}A_{m}$, $A_{n}B_{m}$, $B_{n}B_{m}$:

$\displaystyle\begin{split}A_{n}A_{m}&=\left(\alpha_{n}\psi^{n}+\alpha^{*}_{n}\psi^{-n}\right)\left(\alpha_{m}\psi^{m}+\alpha^{*}_{m}\psi^{-m}\right)\\\
&=\alpha_{n}\alpha_{m}\psi^{n+m}+\alpha^{*}_{n}\alpha_{m}\psi^{m-n}\\\
&\qquad+\alpha_{n}\alpha^{*}_{m}\psi^{n-m}+\alpha^{*}_{n}\alpha^{*}_{m}\psi^{-n-m)}\end{split}$
(70)
$\displaystyle\begin{split}A_{n}B_{m}&=-i\left(\alpha_{n}\psi^{n}+\alpha^{*}_{n}\psi^{-n}\right)\left(\alpha_{m}\psi^{m}-\alpha^{*}_{m}\psi^{-m}\right)\\\
&=-i\left\\{\alpha_{n}\alpha_{m}\psi^{n+m}+\alpha^{*}_{n}\alpha_{m}\psi^{m-n}\right.\\\
&\qquad\left.-\alpha_{n}\alpha^{*}_{m}\psi^{n-m}-\alpha^{*}_{n}\alpha^{*}_{m}\psi^{-n-m)}\right\\}\end{split}$
(71)
$\displaystyle\begin{split}B_{n}B_{m}&=-\left(\alpha_{n}\psi^{n}-\alpha^{*}_{n}\psi^{-n}\right)\left(\alpha_{m}\psi^{m}-\alpha^{*}_{m}\psi^{-m}\right)\\\
&=-\left\\{\alpha_{n}\alpha_{m}\psi^{n+m}-\alpha^{*}_{n}\alpha_{m}\psi^{m-n}\right.\\\
&\qquad\left.-\alpha_{n}\alpha^{*}_{m}\psi^{n-m}+\alpha^{*}_{n}\alpha^{*}_{m}\psi^{-n-m}\right\\}\end{split}$
(72)

Now we have that

$\displaystyle\begin{split}MM_{nm}&=A_{n}A_{m}CC_{nm}+2A_{n}B_{m}CS_{nm}+B_{n}B_{m}SS_{nm}\\\
&=\alpha_{n}\alpha_{m}\widetilde{CC}_{nm}\psi^{n+m}+2\alpha_{n}\alpha^{*}_{m}\widetilde{CS}_{nm}\psi^{n-m}\\\
&\qquad+\alpha^{*}_{n}\alpha^{*}_{m}\widetilde{SS}_{nm}\psi^{-(n+m)}\end{split}$
(73)

where

$\displaystyle\widetilde{CC}_{nm}$
$\displaystyle=(CC_{nm}-SS_{nm})-i\left(CS+CS^{T}\right)_{nm}$ (74)
$\displaystyle\widetilde{CS}_{nm}$ $\displaystyle=(CC_{nm}+SS_{nm})+i\left(CS-
CS^{T}\right)_{nm}$ (75) $\displaystyle\widetilde{SS}_{nm}$
$\displaystyle=(CC_{nm}-SS_{nm})+i\left(CS+CS^{T}\right)_{nm}$ (76)

and for $YM$:

$\displaystyle\begin{split}YM_{k}&=A_{k}YC_{k}+B_{k}YS_{k}\\\
&=\alpha_{k}YC_{k}\psi^{k}+\alpha_{k}^{*}YC_{k}\psi^{-k}\\\
&\qquad-i\left(\alpha_{k}YS_{k}\psi^{k}-\alpha^{*}_{k}YS_{k}\psi^{-k}\right)\\\
&=(YC_{k}-iYS_{k})\alpha_{k}\psi^{k}+(YC_{k}+iYS_{k})\alpha^{*}_{k}\psi^{-k}\\\
&=\alpha_{k}\widetilde{YC}_{k}\psi^{k}+\alpha^{*}_{k}\widetilde{YC}_{k}^{*}\psi^{-k}.\end{split}$
(77)

We also define $YM^{\prime}=\psi^{H}YM$ and $MM^{\prime}=\psi^{2H}MM$, both of
which are polynomials in $\psi$. Their derivatives are

$\displaystyle\partial YM^{\prime}$
$\displaystyle=H\psi^{H-1}YM+\psi^{H}\partial YM$ (78) $\displaystyle\partial
MM^{\prime}$ $\displaystyle=2H\psi^{2H-1}MM+\psi^{2H}\partial MM.$ (79)

A new polynomial condition can then be expressed in terms of $MM^{\prime}$,
$YM^{\prime}$ and their derivatives.

$\displaystyle 0$ $\displaystyle=2MM\partial(YM)-YM\partial(MM)$ (80)
$\displaystyle=\psi^{3H+1}\left(2MM\partial(YM)-YM\partial(MM)\right)$ (81)
$\displaystyle=2\psi^{2H}MM\left(\psi^{H+1}\partial(YM)\right)$
$\displaystyle\qquad-\psi^{H}YM\left(\psi^{2H+1}\partial(MM)\right)$ (82)
$\displaystyle=2MM^{\prime}\left(\psi\partial(YM^{\prime})-H(YM^{\prime})\right)$
$\displaystyle\qquad-
YM^{\prime}\left(\psi\partial(MM^{\prime})-2H(MM^{\prime})\right)$ (83)
$\displaystyle=\psi\left(2MM^{\prime}\partial(YM^{\prime})-YM^{\prime}\partial(MM^{\prime})\right)$
$\displaystyle\qquad-2H\left(MM^{\prime}YM^{\prime}-YM^{\prime}MM^{\prime}\right)$
(84)
$\displaystyle=2MM^{\prime}\partial(YM^{\prime})-YM^{\prime}\partial(MM^{\prime})$
(85)

The last step assumes that $\psi\neq 0$, which is a valid assumption since
$\psi=e^{i\theta_{2}}$ lies on the unit circle for all real $\theta_{2}$.

We solve for the zeros of the polynomial condition defined by Equation 80
using the numpy.polynomial.polyroots function, which solves for the
eigenvalues of the polynomial companion matrix.

Solving for the zeros of a polynomial given a set of coefficients is unstable
in certain cases, since the coefficients are represented as floating point
numbers with finite precision. Thus, we scale the roots by their modulus to
ensure they lie on the unit circle. Alternatively, we could use iterative
schemes such as Newton’s method to improve the estimate of the roots more
robustly, however this requires more computational power and the accuracy of
the roots was not a problem for any of the cases the authors have tested.

###  2.1. Negative amplitude solutions

The model $\hat{y}=\theta_{1}\mathbf{M}_{\theta_{2}}+\theta_{0}$ allows for
$\theta_{1}<0$ solutions. In the original formulation of Lomb-Scargle and in
linear extensions involving multiple harmonics, negative amplitudes translate
to phase differences, since $-\cos{x}=\cos(x-\pi)$ and $-\sin{x}=\sin(x-\pi)$.

However, for non-sinusoidal templates, $\mathbf{M}$, negative amplitudes do
not generally correspond to a phase difference. For example, a detached
eclipsing binary template $\mathbf{M}_{\rm EB}(x)$ cannot be expressed in
terms of a phase-shifted negative eclipsing binary template; i.e.
$\mathbf{M}_{\rm EB}\neq-\mathbf{M}_{\rm EB}(x-\phi)$ for any
$\phi\in[0,2\pi)$.

Negative amplitude solutions found by the fast template periodogram are
usually undesirable, as they may produce false positives for lightcurves that
resemble flipped versions of the desired template, and allowing for
$\theta_{1}<0$ solutions increases the number of effective free parameters of
the model, which lowers the signal to noise, especially for weak signals.

One possible remedy for this problem is to set $P_{\rm FTP}(\omega)=0$ if the
optimal solution for $\theta_{1}$ is negative, but this complicates the
interpretation of $P_{\rm FTP}$. Another possible remedy is, for frequencies
that have a $\theta_{1}<0$ solution, to search for the optimal parameters
while enforcing that $\theta_{1}>0$, e.g. via non-linear optimization, but
this likely will eliminate the computational advantage of FTP over existing
methods.

Thus, we allow for negative amplitude solutions in the model fit and caution
the user to check that the best fit $\theta_{1}$ is positive.

###  2.2. Extending to multi-band observations

####  2.2.1 Multi-phase model

As shown in VanderPlas & Ivezić (2015), the multi-phase periodogram (their
$(N_{\rm base},N_{\rm band})=(0,1)$ periodogram), for any model can be
expressed as a linear combination of single-band periodograms:

$P^{(0,1)}(\omega)=\frac{\sum_{k=1}^{K}\chi^{2}_{0,k}P_{k}(\omega)}{\sum_{k=1}^{K}\chi^{2}_{0,k}}$
(86)

where $K$ denotes the number of bands, $\chi^{2}_{0,k}$ is the weighted sum of
squared residuals between the data in the $k$-th band and its weighted mean
$\left<y\right>$, and $P_{k}(\omega)$ is the periodogram value of the $k$-th
band at the trial frequency $\omega$.

With Equation 86, the template periodogram is readily applicable to multi-band
time series, which is crucial for experiments like LSST, SDSS, Pan-STARRS, and
other current and future photometric surveys.

Other multi-band extensions of the template periodogram are provided in
Appendix A.

###  2.3. Computational requirements

For a given number of harmonics $H$, the task of deriving the polynomial given
in Equation 80 requires $\mathcal{O}(H^{2})$ computations, and finding the
roots of this polynomial requires $\mathcal{O}(H^{3})$ computations. The
degree of the final polynomial is $6H-1$.

When considering $N_{f}$ trial frequencies, the polynomial computation and
root-finding step scales as $\mathcal{O}(H^{3}N_{f})$. The computation of the
sums (Equations 59 – 63) scales as $\mathcal{O}(HN_{f}\log HN_{f})$.

Therefore, the entire template periodogram scales as

$\mathcal{O}(HN_{f}\log HN_{f}+H^{3}N_{f}).$ (87)

However, an important consideration is that the peak width scales inversely
with the number of harmonics $\delta f_{\mathrm{peak}}\propto 1/H$ and so the
number of trial frequencies needed to resolve a peak increases linearly with
the number of harmonics in the template.

The computational scaling factoring in an extra power of $H$ is therefore

$\mathcal{O}(H^{2}N_{f}\log HN_{f}+H^{4}N_{f}).$ (88)

Figure 1.—  Computation time of FTP scaled by $NH$ for different numbers of
harmonics. For $H\lesssim 3$, FTP scales sublinearly in $H$ (possibly due to a
constant overhead per trial frequency, independent of $H$). When $3\lesssim
H\lesssim 11$, FTP scales approximately linearly in $H$, and when $H\gtrsim
11$ FTP approaches the $\mathcal{O}(H^{3})$ scaling limit.

For a fixed number of harmonics $H$, the template periodogram scales as
$\mathcal{O}(N_{f}\log N_{f})$. However, for a constant number of trial
frequencies $N_{f}$, the template algorithm scales as $\mathcal{O}(H^{3})$,
and computational resources alone limit $H$ to reasonably small numbers
$H\lesssim 15$ (see Figure 1).

##  3\. Implementation

Figure 2.—  Template periodograms performed on a simulated eclipsing binary
lightcurve (shown phase-folded in the left-hand plots). The top-most plot uses
only one harmonic, equivalent to a Lomb-Scargle periodogram. Subsequent plots
use an increasing number of harmonics, which produces a narrower and higher
peak height around the correct frequency. For comparison, the multi-harmonic
extension to Lomb-Scargle is plotted in blue, using the same number of
harmonics as the FTP. The Box Least-Squares (Kovács et al., 2002) periodogram
is shown in the final plot.

An open-source implementation of the template periodogram in Python is
available.111https://github.com/PrincetonUniversity/FastTemplatePeriodogram
Polynomial algebra is performed using the numpy.polynomial module (Jones et
al., 2001–). The nfft Python module, 222https://github.com/jakevdp/nfft which
provides a Python implementation of the non-equispaced fast Fourier transform,
is used to compute the necessary sums for a particular time series.

No explicit parallelism is used anywhere in the current implementation,
however certain linear algebra operations in Scipy use OpenMP via calls to
BLAS libraries that have OpenMP enabled.

All timing tests were run on a quad-core 2.6 GHz Intel Core i7 MacBook Pro
laptop (mid-2012 model) with 8GB of 1600 MHz DDR3 memory. The Scipy stack
(version 0.18.1) was compiled with multi-threaded MKL libraries.

###  3.1. Comparison with non-linear optimization

In order to evaluate the accuracy and speed of the template periodogram, we
have included slower alternative solvers within the Python implementation of
the FTP that employ non-linear optimization to find the best fit parameters.

Periodograms computed in Figures 2, 3, and 4 used simulated data. The
simulated data has uniformly random observation times, with Gaussian-random,
homoskedastic, uncorrelated uncertainties. An eclipsing binary template,
generated by fitting a well-sampled, high signal-to-noise eclipsing binary in
the HATNet dataset (BD+56 603) with a 10-harmonic truncated Fourier series.

####  3.1.1 Accuracy

Figure 3.—  Comparing accuracy between previous methods that rely on non-
linear optimization at each trial frequency with the fast template periodogram
described in this paper. Both methods are applied to the same simulated data
as shown in Figure 2. The FTP consistently finds more optimal template fits
than those found with non-linear optimization, which do not guarantee
convergence to a globally optimal solution. The FTP solves for the optimal fit
parameters directly, and therefore is able to achieve greater accuracy than
template fits done via non-linear optimization. Figure 4.—  Comparing the
template periodogram calculated with $H=10$ harmonics to the template
periodogram using a smaller number of harmonics $H<10$. The template and data
used to perform the periodogram calculations are the same as those shown in
Figure 2.

For weak signals or signals folded at the incorrect trial period, there may be
a large number of local $\chi^{2}$ minima in the parameter space, and thus
non-linear optimization algorithms may have trouble finding the global
minimum. The FTP, on the other hand, solves for the optimal parameters
directly, and thus is able to recover optimal solutions even when the signal
is weak or not present.

Figure 3 illustrates the accuracy improvement with FTP. Many solutions found
via non-linear optimization are significantly suboptimal compared to the
solutions found by the FTP.

Figure 4 compares FTP results obtained using the full template $(H=10)$ with
those obtained using smaller numbers of harmonics. The left-most plot compares
the $H=1$ case (weighted Lomb-Scargle), which, as also demonstrated in Figure
2, illustrates the advantage of the template periodogram for known, non-
sinusoidal signal shapes.

####  3.1.2 Computation time

Figure 5.—  Computation time of FTP compared with alternative techniques that
use non-linear optimization at each trial frequency. _Left_ : timing for the
case when $N_{f}=~{}12N_{\rm obs}$, i.e. the cadence of the observations is
constant. _Right_ : timing for the case when $N_{f}$ is fixed, i.e. the
baseline of the observations are constant. Non-linear optimization techniques
scale as $\mathcal{O}(N_{f}N_{\rm obs})$ while the FTP scales as
$\mathcal{O}(HN_{f}\log HN_{f}+N_{f}H^{3})$, where $H$ is the number of
harmonics needed to approximate the template.

FTP scales asymptotically as $\mathcal{O}(N_{f}H\log N_{f}H)$ with respect to
the number of trial frequencies, $N_{f}$ and as $\mathcal{O}(N_{f}H^{3})$ with
respect to the number of harmonics in which the template is expanded, $H$. For
a given resolving power, however, there is an additional factor of $H$ in each
of these terms due to the number of trial frequencies necessary to resolve a
periodogram peak being proportional to $H$. However, for reasonable cases
($N_{f}\lesssim 10^{120}$ when $H=5$) the computation time is dominated by
computing polynomial coefficients and root finding, both of which scale
linearly in $N_{f}$.

The number of trial frequencies needed for finding astrophysical signals in a
typical photometric time series is

$N_{f}=1.75\times 10^{6}\left(\frac{H}{1}\right)\left(\frac{{\rm
baseline}}{10~{}{\rm
yrs}}\right)\left(\frac{\alpha}{5}\right)\left(\frac{15~{}{\rm mins}}{P_{\rm
min}}\right)$ (89)

where $\alpha$ represents the “oversampling factor,” $\Delta f_{\rm
peak}/\Delta f$, where $\Delta f_{\rm peak}\sim 1/{\rm baseline}$ is the
typical width of a peak in the periodogram and $\Delta f$ is the frequency
spacing of the periodogram.

Extrapolating from the timing of a test case (500 observations, 5 harmonics,
15,000 trial frequencies), the summations account for approximately 5% of the
computation time when $N_{f}\sim 10^{6}$. If polynomial computations and root-
finding can be improved to the point where they no longer dominate the
computation time, this would provide an order of magnitude speedup over the
current implementation.

Figure 5 compares the timing of the FTP with that of previous methods that
employ non-linear optimization. For the case when $N_{f}\propto N_{\rm obs}$,
FTP achieves a factor of 3 speedup for even the smallest test case (15
datapoints), while for larger cases ($N\sim 10^{4}$) FTP offers 2-3 orders of
magnitude speed improvement. For the constant baseline case, FTP is a factor
of $\sim 2$ faster for the smallest test case and a factor of $\sim 20$ faster
for $N_{\rm obs}\sim 10^{4}$. Future improvements to the FTP implementation
could further improve speedups by 1-2 orders of magnitude over non-linear
optimization.

##  4\. Discussion

Template fitting is a powerful technique for accurately recovering the period
and amplitude of objects with _a priori_ known lightcurve shapes. It has been
used in the literature by, e.g. Stringer et al. (2019); Sesar et al. (2016,
2010), to analyze RR Lyrae in the SDSS, PS1, and DES datasets, where it has
been shown to produce purer samples of RR Lyrae at a given completeness. The
computational cost of current template fitting algorithms, however, limits
their application to larger datasets or with a larger number of templates.

We have presented a novel template fitting algorithm that extends the Lomb-
Scargle periodogram (Lomb, 1976; Scargle, 1982; Barning, 1963; Vaníček, 1971)
to handle non-sinusoidal signals that can be expressed in terms of a truncated
Fourier series with a reasonably small number of harmonics ($H\lesssim 10$).

The fast template periodogram (FTP) asymptotically scales as
$\mathcal{O}(N_{f}H^{2}\log N_{f}H^{2}+N_{f}H^{4})$, while previous template
fitting algorithms such as the one used in the gatspy library (VanderPlas,
2016), scale as $\mathcal{O}(N_{f}N_{\rm obs})$. However, the FTP effectively
scales as $\mathcal{O}(N_{f}H^{4})$, since the time needed to compute
polynomial coefficients and perform zero-finding dominates the computational
time for all practical cases ($N_{f}\lesssim 10^{120}$). The $H^{4}$ scaling
effectively restricts templates to those that are sufficiently smooth to be
explained by a small number of Fourier terms.

FTP also improves the accuracy of previous template fitting algorithms, which
rely on non-linear optimization at each trial frequency to minimize the
$\chi^{2}$ of the template fit. The FTP routinely finds superior fits over
non-linear optimization methods.

An open-source Python implementation of the FTP is available at
GitHub.333https://github.com/PrincetonUniversity/FastTemplatePeriodogram The
current implementation could likely be improved by:

  1. 1.

Improving the speed of the polynomial coefficient calculations and the zero-
finding steps. This could potentially yield a speedup of $\sim 1-2$ orders of
magnitude over the current implementation.

  2. 2.

Exploiting the embarassingly parallel nature of the FTP using GPU’s.

For a constant baseline, the current implementation improves existing methods
by factors of a $\sim$few for lightcurves with $\mathcal{O}(100)$
observations, and by an order of magnitude or more for objects with more than
1,000 observations. These improvements, taken at face value, are not enough to
make template fitting feasible on LSST-sized datasets. However, optimizing the
polynomial computations could yield a factor of $\sim 25-100$ speedup over the
current implementation, which would make the FTP 1-3 orders of magnitude
faster than alternative techniques.

##  Appendix A Shared-phase multi-band template periodogram

We derive a multi-band extension for the template periodogram for data taken
in $K$ filters, with $N_{k}$ observations in the $k$-th filter. We use the
same model as the one described in Sesar et al. (2016) in order to illustrate
the applicability of the template periodogram to more sophisticated scenarios.

The $i$-th observation in the $k$-th filter is denoted $y^{(k)}_{i}$. We wish
to fit a periodic, multi-band template
$\mathbf{M}=(\mathbf{M}^{(1)},\mathbf{M}^{(2)},...,\mathbf{M}^{(K)}):[0,2\pi)\rightarrow\mathbb{R}^{K}$
to all observations. We assume the same model used by Sesar et al. (2016),
which assumes the relative amplitudes, phase shifts, and offsets are shared
across bands:

$\hat{y}^{(k)}(t|\theta)=\theta_{1}\mathbf{M}^{(k)}(\omega
t-\theta_{2})+\theta_{3}+\lambda^{(k)}$ (A1)

where $\lambda^{(k)}$ is a fixed relative offset for band $k$. The $\chi^{2}$
for this model is

$\chi^{2}=\sum_{k=1}^{K}\sum_{i=1}^{N_{k}}w^{(k)}_{i}\left(y^{(k)}_{i}-\hat{y}^{(k)}_{i}\right)^{2}$
(A2)

To make things simpler, we can set the $\lambda^{(k)}$ values to 0 simply by
subtracting them off from our observations; this means we take all
$y^{(k)}_{i}\rightarrow y^{(k)}_{i}-\lambda^{(k)}$. We have that

$\chi^{2}=\sum_{k=1}^{K}W^{(k)}\chi^{2}_{k}$ (A3)

Where $W^{(k)}\equiv\sum_{i=1}^{N_{k}}w^{(k)}_{i}$ and
$\sum_{k=1}^{K}W^{(k)}=1$. This means that the system of equations reduces to:

$0=\frac{\partial\chi^{2}}{\partial\theta_{1}}=\sum_{k=1}^{K}W^{(k)}\left(\widehat{YM}^{(k)}-\theta_{1}\widehat{MM}^{(k)}-\theta_{3}\overline{M}^{(k)}\right)$
(A4)

for the $\theta_{1}$ parameter, where $\widehat{YM}^{(k)}$,
$\widehat{MM}^{(k)}$, $\overline{M}^{(k)}$ are values from the single band
case computed for each band individually, holding $W^{(k)}=1$ for each band.
That is:

$\displaystyle\widehat{YM}^{(k)}$
$\displaystyle\equiv\frac{1}{W^{(k)}}\sum_{i=1}^{N_{k}}w^{(k)}_{i}y^{(k)}_{i}\mathbf{M}_{\theta_{2}}^{(k)}(\omega
t_{i})$ (A5) $\displaystyle\widehat{MM}^{(k)}$
$\displaystyle\equiv\frac{1}{W^{(k)}}\sum_{i=1}^{N_{k}}w^{(k)}_{i}\left(\mathbf{M}_{\theta_{2}}^{(k)}(\omega
t_{i})\right)^{2}$ (A6) $\displaystyle\overline{M}^{(k)}$
$\displaystyle\equiv\frac{1}{W^{(k)}}\sum_{i=1}^{N_{k}}w^{(k)}_{i}\mathbf{M}_{\theta_{2}}^{(k)}(\omega
t_{i})$ (A7)

For the offset $\theta_{3}$, we have

$0=\frac{\partial\chi^{2}}{\partial\theta_{3}}=\sum_{k=1}^{K}W^{(k)}\left(\bar{y}^{(k)}-\theta_{1}\overline{M}^{(k)}-\theta_{3}\right)$
(A8)

Where $\bar{y}^{(k)}$ is the weighted-mean for the $k$-th band
($\bar{y}^{(k)}\equiv(1/W^{(k)})\sum_{i=1}^{N_{k}}w^{(k)}_{i}y^{(k)}_{i}$),
again with $y^{(k)}_{i}\rightarrow y^{(k)}_{i}-\lambda^{(k)}$.

So if we redefine the quantities $\widehat{YM}$, $\widehat{MM}$,
$\overline{M}$, and $\bar{y}$ as weighted averages across the bands, i.e.
$\widehat{YM}=\sum_{k=1}^{K}W^{(k)}\widehat{YM}^{(k)}$, etc., the solution for
the optimal parameters has the same form:

$\displaystyle\theta_{1}$ $\displaystyle=\left(\frac{YM}{MM}\right)$ (A9)
$\displaystyle\theta_{3}$
$\displaystyle=\bar{y}-\overline{M}\left(\frac{YM}{MM}\right)$ (A10)

which means that the model for the $k$-th band is

$\hat{y}^{(k)}=\bar{y}+\left(\frac{YM}{MM}\right)\left(M^{(k)}_{i}-\overline{M}\right)$
(A11)

The form of the periodogram has the same form as the single-band case:

$\displaystyle\chi^{2}$
$\displaystyle=\sum_{k=1}^{K}\sum_{i=1}^{N_{k}}w^{(k)}_{i}\left(y^{(k)}_{i}-\hat{y}^{(k)}_{i}\right)^{2}$
(A12)
$\displaystyle=\sum_{k=1}^{K}\sum_{i=1}^{N_{k}}w^{(k)}_{i}\left((y^{(k)}_{i}-\bar{y})-\left(\frac{YM}{MM}\right)(M_{i}^{(k)}-\overline{M})\right)^{2}$
(A13)
$\displaystyle=\sum_{k=1}^{K}\sum_{i=1}^{N_{k}}w^{(k)}_{i}\left((y^{(k)}_{i}-\bar{y})^{2}-2\left(\frac{YM}{MM}\right)(M_{i}^{(k)}-\overline{M})(y^{(k)}_{i}-\bar{y})+\left(\frac{YM}{MM}\right)^{2}(M_{i}^{(k)}-\overline{M})^{2}\right)^{2}$
(A14) $\displaystyle=YY-2\frac{(YM)^{2}}{MM}+\frac{(YM)^{2}}{MM}$ (A15)
$\displaystyle=YY-\frac{(YM)^{2}}{MM}.$ (A16)

Since there is a single shared offset between the bands (i.e. we assume the
mean magnitude is the same in all bands after subtracting $\lambda^{(k)}$),
the variance for the signal, $YY$, is not
$\sum_{k=1}^{K}\sum_{i=1}^{N_{k}}w^{(k)}_{i}(y^{(k)}_{i}-\bar{y}^{(k)})^{2}$
but $\sum_{k=1}^{K}\sum_{i=1}^{N_{k}}w^{(k)}_{i}(y^{(k)}_{i}-\bar{y})^{2}$.

Construction of the polynomial for the multi-band case can be performed by
first computing the polynomial expression for
$\psi^{2H}\widehat{MM}=\sum_{k=1}^{K}\sum_{i=1}^{N_{k}}w^{(k)}_{i}\left(\psi^{H}M^{(k)}_{i}\right)^{2}$
and a separate polynomial expression for $\psi^{H}\overline{M}$, which can
then be squared and subtracted from $\psi^{2H}\widehat{MM}$ to find
$MM^{\prime}=\psi^{2H}MM$.

For $YM^{\prime}=\psi^{H}YM$, we merely compute $\psi^{H}\widehat{YM}^{(k)}$
for each band, take the weighted average of the polynomial coefficients for
all bands
($\psi^{H}\widehat{YM}=\sum_{k}^{K}W^{(k)}\psi^{H}\widehat{YM}^{(k)}$) and
subtract the $\psi^{H}\bar{y}\overline{M}$ polynomial to get $YM^{\prime}$.

After computing the polynomials $MM^{\prime}$ and $YM^{\prime}$, the
polynomial in Equation 80 can be computed quickly and the following steps for
finding the optimal model parameters is the same as in the single band case.

The computational complexity of this model scales as $\mathcal{O}(KN_{f}H\log
N_{f}H+N_{f}H^{3})$ where $K$ is the number of filters. Since the polynomial
zero-finding step is the limiting computation in most real-world applications,
the multi-band template periodogram corresponding to the Sesar et al. (2016)
model should not be significantly more computationally intensive than the
single-band case.

Joel Hartman and GB acknowledge support from NASA grant NNX17AB61G. JTV is
supported by the University of Washington eScience Institute, with funding
from the Alfred P. Sloan Foundation, the Gordon and Betty Moore Foundation,
and the Washington Research Foundation.

## References

  * Bakos et al. (2004) Bakos, G., Noyes, R. W., Kovács, G., et al. 2004, PASP, 116, 266
  * Barning (1963) Barning, F. J. M. 1963, Bull. Astron. Inst. Netherlands, 17, 22
  * Bretthorst et al. (1988) Bretthorst, G. L., Hung, C.-C., D’Avignon, D. A., & Ackerman, J. J. H. 1988, Journal of Magnetic Resonance, 79, 369
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, ArXiv e-prints, arXiv:1612.05560
  * Cincotta et al. (1995) Cincotta, P. M., Mendez, M., & Nunez, J. A. 1995, ApJ, 449, 231
  * Cooley & Tukey (1965) Cooley, J. W., & Tukey, J. W. 1965, Math. Comput., 19, 297
  * Dutt & Rokhlin (1993) Dutt, A., & Rokhlin, V. 1993, SIAM J. Sci. Comput., 14, 1368
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Graham et al. (2013a) Graham, M. J., Drake, A. J., Djorgovski, S. G., Mahabal, A. A., & Donalek, C. 2013a, MNRAS, 434, 2629
  * Graham et al. (2013b) Graham, M. J., Drake, A. J., Djorgovski, S. G., et al. 2013b, MNRAS, 434, 3423
  * Hernitschek et al. (2016) Hernitschek, N., Schlafly, E. F., Sesar, B., et al. 2016, ApJ, 817, 73
  * Huijse et al. (2012) Huijse, P., Estevez, P. A., Protopapas, P., Zegers, P., & Principe, J. C. 2012, IEEE Transactions on Signal Processing, 60, 5135
  * Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python, [Online; accessed 2017-01-19]
  * Keiner et al. (2009) Keiner, J., Kunis, S., & Potts, D. 2009, ACM Trans. Math. Softw., 36, 19:1
  * Kovács et al. (2002) Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369
  * Leroy (2012) Leroy, B. 2012, A&A, 545, A50
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, ArXiv e-prints, arXiv:0912.0201
  * Palmer (2009) Palmer, D. M. 2009, ApJ, 695, 496
  * Press & Rybicki (1989) Press, W. H., & Rybicki, G. B. 1989, ApJ, 338, 277
  * Rasmussen & Williams (2005) Rasmussen, C. E., & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) (The MIT Press)
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
  * Schwarzenberg-Czerny (1996) Schwarzenberg-Czerny, A. 1996, ApJ, 460, L107
  * Sesar et al. (2010) Sesar, B., Ivezić, Ž., Grammer, S. H., et al. 2010, ApJ, 708, 717
  * Sesar et al. (2016) Sesar, B., Hernitschek, N., Mitrović, S., et al. 2016, ArXiv e-prints, arXiv:1611.08596
  * Stringer et al. (2019) Stringer, K. M., Long, J. P., Macri, L. M., et al. 2019, AJ, 158, 16
  * VanderPlas (2016) VanderPlas, J. 2016, gatspy: General tools for Astronomical Time Series in Python, Astrophysics Source Code Library, ascl:1610.007
  * VanderPlas (2018) VanderPlas, J. T. 2018, ApJS, 236, 16
  * VanderPlas & Ivezić (2015) VanderPlas, J. T., & Ivezić, Ž. 2015, ApJ, 812, 18
  * Vaníček (1971) Vaníček, P. 1971, Ap&SS, 12, 10
  * Zechmeister & Kürster (2009) Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577


# Evaluating Robot Posture Control and Balance by Comparison to Human Subjects
using Human Likeness Measures

Vittorio Lippi1 , Christoph Maurer1 , Thomas Mergner1
1University Clinic of Freiburg, Neurology, Freiburg, Germany
{vittorio.lippi , christoph.maurer<EMAIL_ADDRESS>
https://orcid.org/0000-0001-5520-8974 https://orcid.org/0000-0001-9050-279X
https://orcid.org/0000-0001-7231-164X

###### Abstract

Posture control and balance are basic requirements for a humanoid robot
performing motor tasks like walking and interacting with the environment. For
this reason, posture control is one of the elements taken into account when
evaluating the performance of humanoids. In this work, we describe and analyze
a performance indicator based on the comparison between the body sway of a
robot standing on a moving surface and the one of healthy subjects performing
the same experiment. This approach is here oriented to the evaluation of human
likeness. The measure is tested with three human-inspired humanoid posture
control systems, the independent channel (IC), the disturbance identification
and compensation (DEC), and the eigenmovement (EM) control. The potential and
the limitations connected with such human-inspired humanoid control mechanisms
are then discussed.

##  1 INTRODUCTION

The benchmarking of humanoid performance is gaining interest in the research
community [Conti et al., 2018, Torricelli et al., 2014, Torricelli et al.,
2020, Lippi et al., 2020c, Lippi et al., 2019]. The performance of a humanoid
is a complex subject that covers several issues, e.g. sensor fusion, cognitive
and motor aspects, mechanics, and energy efficiency. In particular, a recent
European project, EUROBENCH, is proposing to implement standard and repeatable
experimental procedures to evaluate and compare the performance of different
robots. This work describes one specific posture control performance indicator
implemented within the project, evaluating human-likeness based on the
similarity between responses to external disturbances, more are described in
[Lippi et al., 2019] and [Lippi et al., 2020c]. Balance has been taken into
account because falling is one of the typical reasons of failure for humanoids
[Atkeson et al., 2015, Atkeson et al., 2018, Guizzo and Ackerman, 2015]. In
particular, the aspect of human likeliness is studied in this work. A formal
and unanimous definition of human likeliness is still missing, although the
concept is relevant both for robotics and neuroscience [Torricelli et al.,
2014].

Figure 1: Posture control experiment with a moving support surface, a 6-Dof
Stuart Platform. On the left, a human subject with active markers visible on
the back (blue triangular plates) is shown. On the right, a robot standing on
the same platform.

Human-like behavior is expected to entail some desirable features such as
mechanical compliance (that is important for safety) and low energy
consumption. Both of these features are reasonably associated with the low
feedback loop gain imposed on biological systems in face of the presence of
neural delay. In robotics the delay associated with the transport of signals
is negligible, but the robustness in face of delays may be important in
limiting the required computational resources needed to implement the real-
time control [Ott et al., 2016]. Although being often considered important,
human-likeness at the state of the art is defined in terms of perception from
the point of view of human observers [von Zitzewitz et al., 2013, Oztop et
al., 2005, Abel et al., 2020] or identified with the presence of a specific
feature, e.g. reproducing human trajectories [Kim et al., 2009] or exhibiting
compliance [Colasanto et al., 2015]. The measure proposed here is based on
body kinematics to make the evaluation of human likeness repeatable and
objective instead of designed ad hoc and subjective. Specifically, the
evaluation is based on a data-set of results from human experiments. The
experimental data consist of the body sway of human body segments induced by
an external stimulus, here specifically the tilt of the support surface in the
body sagittal plane. The measured body sway is characterized as a frequency
response function, FRF, i.e. an empirical transfer function between the PRTS
stimulus and the response computed on specific frequency points. The data-set
has been produced with healthy subjects. The comparison aims to assess the
similarity in the balancing behavior between the robot and an average healthy
subject. Such measure does not make explicit assumptions about specific
advantages of the human behavior, in contrast to the performance indexes
proposed in other works that identify some goal like minimum torque or energy
consumption in robots [Kashiri et al., 2016], or a specific problem in human
subjects [Singh, 2016]. The implications of the proposed measure and the
limitations connected with human likeliness will be discussed in the
conclusions.

##  2 MATERIALS AND METHODS

Comparing body sway profiles FRFs has been chosen as a base for the
benchmarking because such analysis has been studied in several publications
with human subjects. This provides a reference for comparison and tools for
analysis. There are several reasons why the body sway induced by support
surface tilt has been repeatedly used to study human posture control. As a
repeatable stimulus it can be used to formally characterize the behavior in
terms of responses to a well defined input. Furthermore, the tilting support
surface requires the balancing mechanism to integrate vestibular input and
proprioception (and vision, with eyes open), and hence it is well apt to study
human sensor fusion.

Figure 2: PRTS spectrum and averaging process. The magnitude of the discrete
Fourier transform of the PRTS is shown in A. An empirical transfer function
from a trial is shown in B, and the associated 11-component FRF are shown in
C. Note that, as the FRF is averaged in the complex domain the magnitude of
the average shown in the plot is not the average of the magnitudes. The bands
on the background show the frequency ranges over which the spectrum is
averaged: white and dark grey represent ranges associated with groups of
frequencies. The sets of frequencies are overlapping, with light gray bands
belonging to both the two contiguous groups, and a sample on the transition
between two bands belongs to both the groups. The weights in D represents the
values of $P(f_{x})$ in eq. 4.

The sample data-set here considered consists of data from human posture
control experiments. A group of 38 young healthy subjects serves as a
reference for the human-likeness criterion. The raw data include information
like age, body weight, height, and the height of the different markers. In the
experiments, the subjects were presented with a stimulus consisting of a tilt
of the support surface in the body sagittal plane, while the recorded output
was the body sway. The typical set-up is shown in Fig. 1, the tracking was
performed using active markers (Optotrak 3020; Waterloo, ON, Canada), attached
to subjects’ hips and shoulders and to the platform. A PC with custom-made
programs was used to generate the support surface tilt in the body-sagittal
plane. The marker positions were recorded at 100 Hz using software written in
LabView (National Instruments; Austin, TX, United States). The profile used
for the stimulus is a pseudorandom ternary signal, PRTS [Peterka, 2002]. The
peak-to-peak amplitude was set to $1^{\circ}$. The amplitude is rather small
compared to what a healthy subject can withstand. Usually in similar studies
the tilt may go up to 8 degrees and more [Peterka, 2002, Assländer et al.,
2015]. This was motivated by the aim to provide a data-set that could be
compared safely with elderly subjects and patients and, in the specific case
of robotic benchmarking, can be used to characterize the behavior of the robot
without the risk of making the robot fall (and potentially break).
The performance indicator is here a measure of similarity with human behavior.
The body sway profiles, i.e. the angular sway of the body COM with respect to
the ankle joint, is used to characterize and compare the responses. The
comparison is defined in terms of the norm of the difference between frequency
response functions (FRFs). The PRTS power-spectrum has a profile with peaks at
$f_{peak}$ separated by ranges of frequencies with no power [Jilk et al.,
2014, Peterka, 2002, Lippi et al., 2020b]. The response is considered for a
specific set of frequencies where the PRTS spectrum has peaks:
$\mathbf{f_{peak}}$=[0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85,
0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85, 1.95, 2.05, 2.15,
2.25, 2.35, 2.45] Hz. Such a discrete spectrum is then transformed into a
vector of 11 components by averaging the FRF over neighboring frequencies as
illustrated in Fig. 2:

$f_{x(k)}=\frac{\sum_{i\in B_{k}}f_{peak(i)}}{N_{k}}$ (1)

where $k$ and $i$ are the indexes of the components of the frequency vectors,
$B_{k}$ is a set of $N_{k}$ frequencies averaged to obtain the $k^{th}$
sample. The $B_{k}$ are shown in Fig. 2 as white and gray bands (notice that
the bands are overlapping). Similarly the Fourier transform of the PRTS
$P(f_{x})$ and the Fourier transform of the responses are averaged over the
bands $B_{k}$ before computing the FRFs. The final representation of the FRF
is a function of the 11 frequencies $f_{x}$ =[0.1, 0.3, 0.6, 0.8, 1.1, 1.4,
1.8, 2.2, 2.7, 3.5, 4.4]111In a previous description of the system in [Lippi
et al., 2020c] the frequency vector had 16 components as proposed in other
works using PRTS, e.g. [Peterka, 2002, Goodworth and Peterka, 2018, Lippi et
al., 2020b, Lippi et al., 2020a]. In this work, we decided for a shorter
version of the signal which is considered easier for human subjects (sometimes
with elderly patients) and convenient for robotics experimenters. The discrete
Fourier transform of the signal was consequently shorter, as well as the
resulting FRFs.. The choice of the frequencies in $B_{k}$ and their
overlapping follows the method described in [Peterka, 2002] but here is
adapted to a vector of 11 frequencies. The rationale behind such a choice was
to get a representation with the frequencies equally spaced on logarithmic
scale, which is often used in posturography papers to present the FRF, with
the overlapping providing a smoothing effect. In detail, the FRF is computed
from the 11 components of the Fourier transform of the input $U$, and the
output $Y$ as

$FRF=\frac{G_{UY}}{G_{U}}$ (2)

where $G_{UY}=U^{*}\odot Y$ and $G_{U}=U^{*}\odot U$ are empirical estimations
of the cross power spectrum and the input power spectrum (“$\odot$” is the
Hadamard, element-wise product). The peaks of the PRTS power-spectrum have
larger values at lower frequencies [Jilk et al., 2014]. This implies a better
signal-to-noise ratio for the first components. A weighting vector
$\mathbf{w}$ based on $P(f_{x})$ then defined in a similar way to eq. 1, but
considering the power

$w_{k}=\sqrt{\sum_{i\in B_{k}}||P(f_{peak(i)})||^{2}}$ (3)

The distance between two FRFs is defined and the norm of the difference
weighted by the precision matrix, i.e. the inverse of the covariance matrix
$\Sigma$, computed for the data-set of normal subjects. Before doing this the
FRF is expanded into a vector with the real and imaginary components as
separated elements, i.e. 22 components. This together with the foretold
weighting leads to the definition of the norm:

$D=\sqrt{S\Delta^{T}\Sigma^{-1}\Delta S}$ (4)

where $S=diag([\mathbf{w},\mathbf{w}])$ is the diagonal matrix representing
the re-weighting due to the power-spectrum, doubled to cover the 22 elements,
and $\Delta$ is the difference between the two FRFs expanded to 22 components.
This approach does not require model identification because it is performed on
the basis of the data. The comparison can be performed between the tested
robot and the average of the groups of humans (healthy or with special
deficient conditions) or between two single samples to quantify how much two
robots differ from each other. The score of human-likeness is obtained
comparing the sample with the mean of the human sample set $\mathbf{\mu}$ so
that $\Delta=FRF-\mathbf{\mu}$. The parameters $\mathbf{\mu}$ and $\Sigma$
defining the score are shown in Fig. 3.

Figure 3: Parameters $\Sigma$ (top) and $\mathbf{\mu}$ (below) defining the
score. Notice how the values on the diagonal corresponding to the highest
frequencies are associated with high values (less variance in the dataset).
This would make the metrics particularly sensitive to accidental changes in
such components (e.g. due to noise). This disadvantage is compensated by the
weighting profiles in Fig. 2 D, which associate almost zero weight to high-
frequency components. Figure 4: Overview of the tested three control systems.
The notation is slightly different to match the one proposed by the respective
authors in previous works. \tinyA⃝ The single inverted pendulum mode for the
IC, characterized by the body orientation in space $\alpha_{BS}$ and the angle
between body and foot $\alpha_{BF}$. \tinyB⃝. The IC control scheme as
presented in [Pasma et al., 2018]. The three weights $W_{prop}$, $W_{vest}$,
and $W_{vis}$ are associated with proprioception ($\alpha_{BF}$), vestibular
signal ($\alpha_{BS}$) and vision ($\alpha_{BS}$) respectively. The lumped
delay $e^{\Delta t\cdot s}$ represents all the motor and sensory delays in the
loop. The positive forcre feedback includes a gain and a first-order low pass
filter: $K_{F}/(T_{f}s+1)$. It is included to explain the observed behavior at
low frequencies [Peterka, 2003]. \tinyC⃝ The DEC, implemented as double
inverted pendulum (DIP) [Hettich et al., 2014, Hettich et al., 2013, Hettich
et al., 2015], with two control modules: the “Hip Module” controlling the
orientation of the trunk in space $\alpha_{TS}$, and the “Ankle Module”
controlling the orientation of the body CoM in space $\alpha_{BS}$. \tinyD⃝
The scheme of a module of the DEC control. The system includes a passive
control loop (PD) and a neural controller (PID) implementing the servo loop
and the compensation of estimated disturbances as a störgrößenaufschaltung,
i.e. feedforward compensation based on sensory input. There is a lumped delay
$e^{\Delta t\cdot s}$ in each of the modules. The vestibular input vest* for
the ankle module uses a sensor fusion derived information (vestibular +
proprioceptive) to reconstruct the orientation of legs in space. The
vestibular input provides the orientation of the head in space (and hence of
trunk in space) directly. \tinyE⃝ The decomposition of an arbitrary movement
($q_{1}$ = ankle angle, $q_{2}$ = hip angle) in two EMs.\tinyF⃝ The EM control
scheme.

Three control systems are tested: the independent channel model (IC), the
disturbance estimation and compensation (DEC), and the eigenmovements concept
(EM). Next we provide a brief description of the concepts.
1) The IC model is a simplified descriptive linear model for human upright
posture control. The model consists of a single inverted pendulum controlled
by a feedback mechanism with a PD controller and a time delay. The sensory
integration mechanism consists of a weighted sum of the sensory contributions.
The weights are constrained to unity [Peterka, 2002]. The three channels are:
angle between body and foot perceived through proprioception, body-in-space
orientation perceived by the vestibular system, and body-in-space orientation
perceived by vision. The time delay represents all the delays in the control
loop. In order to simulate a robot control experiment similar to the ones
performed with the robot Posturob II in [Pasma et al., 2018] the system has
been simulated in “closed-eyes” conditions. The coefficient $W_{vis}$ is thus
set to zero.
2) The DEC model was proposed as a model of human posture control in steady
state [Mergner et al., 2003]. It was implemented into a humanoid control
system [Lippi and Mergner, 2017, Lippi et al., 2016] and tested on robotic
platforms [Hettich et al., 2014, Lippi, 2018, Ott et al., 2016, Zebenay et
al., 2015]. The DEC exploits sensor-fusion-derived internal reconstructions of
the external disturbances having impact on body posture and balance. The
reconstructed disturbances are compensated using a servo controller. The model
is nonlinear because it includes sensory thresholds and the disturbance
estimators are non-linear functions. The control loop includes a lumped delay
(one for each joint). Its formulation leads to a modular control system with
one control module for each degree of freedom as shown in Fig. 4 (C).
3) The EM [Alexandrov et al., 2005] as implemented in [Alexandrov et al.,
2017]. This control system of bio-mechanical inspiration is based on the
diagonalization of the inertial matrix $B_{0}$ describing the linearized body
dynamics, i.e.

$B_{0}\ddot{\mathbf{q}}-G_{0}\mathbf{q}=\tau_{con}$ (5)

where $\mathbf{q}$ is the vector of joint angles, $G_{0}$ is the gravity
matrix, and $\tau_{con}$ the control torques applied to the joints. EMs have
independent dynamics, i.e. PD controller is implemented for each EM. In the
original formulation the body kinematics is expressed in the joint angles
$q_{1}$ and $q_{2}$ of Fig. 4 E. The control of balance in space is here
implemented by integrating a vestibular system (technically an IMU, inertial
measurement unit) and defining $q_{1}$ as leg-in-space sway angle (i.e. the
sum of the hip angle and vestibular trunk in space orientation). When
simulating the human system, there are different neural delays for the ankle
and trunk control loops. The additional delays in Fig. 4 F are added to make
the delays the same for both the joints, whic here is assumed as a condition
for decoupling the dynamics of the EM.

Figure 5: The hip joint of the robot Posturob II blocked by a screw (left) and
the trunk to leg sway over time in simulation in degrees (right).

The parameters for the simulations are based on the experiments previously
performed with the robot Posturob II using an IC, DEC, and EM approach as
presented in [Pasma et al., 2018], [Hettich et al., 2014], and [Alexandrov et
al., 2017], respectively. The experiments are simulated in MATLAB/Simulink®.
The robot is simulated using the double inverted model from [AlBakri, 2008].
In order to implement the SIP control for the IC, the hip stiffness is set to
a large value, $\approx 10^{5}$ Nm, leading to virtually no movement ($\leq
10^{-3}$ deg, as shown in Fig. 5). This is similar to what was implemented in
[Pasma et al., 2018] blocking the hip joint of the robot with a screw. This
choice of simulation parameters aims to reproduce the robot behavior. The
parameters of the models can be optimized to reproduce the human FRF our in
simulation as previously done in posture control studies, e.g. [Assländer et
al., 2015]. While this is useful to discuss the parameters themselves, e.g. to
estimate the relative weight of different sensory channels in the control loop
in different conditions, it is not obviously leading to a human-like behavior
when applied on real robots. The responses are compared with eyes-closed human
experiments and hence the visual signal is not considered: a trunk-in-space
orientation is assumed to be available in the simulations with no specific
modeling of the dynamics of the vestibular system (which is, in general,
faster but less accurate than the visual signal). This is because, at the
state of the art, humanoids are not (yet) using visual input for balance
control. This may change in the future considering the current advancements in
computer vision and their integration in humanoid robotics. Data-sets for the
eyes-open condition are however available for future benchmarking.

##  3 RESULTS

The simulation with the three control systems produced the FRFs shown in Fig.
6 together with the average of the human FRFs in the reference data-set. The
results obtained were:

control | score | CDF
---|---|---
IC | 2.2034 | 75.2212 %
DEC | 2.7318 | 87.6106 %
EM | 3.9616 | 97.3451 %

The values in the third column give estimates of the cumulative distribution
function (CDF) of the data-set at the given score. It is computed by counting
the fraction of samples with a score smaller than the one produced in the
specific simulation. The comparison in Fig. 6 shows that all three robotic
controllers exhibit a smaller gain than the one observed on average in humans
(especially the EM) and differences of phase. In Fig. 7 the result for the
three controllers is compared with the distribution of the human data-set.

##  4 DISCUSSION, CONCLUSIONS AND FUTURE WORK

The main contribution of the present work is a set of data from human
experiments that can be used as a benchmark for human-like posture control to
be compared with robots, and the evaluation of a metric proposed for such
comparison.
The tests for the bio-inspired control systems showed that the three robot
control systems produced small gains compared to humans: according to the
proposed measure they were “less human-like” than most of the samples in the
data-set. Notice that, as explained in the methods section, the parameters of
the models can be optimized to create a FRF that is more similar to the human
average (as it is actually done in several posture control studies that
include a control model). Here, the aim of the simulation was not to evaluate
the three models in general (the conclusion here is not that the IC is more
human-like from an absolute point of view) or to have a “fair competition”
between the three control concepts: that would require to perform the test on
the real robot or, at least, optimize the control parameters based on the
controlled system. The main idea was to get real models as they were used in
real robot experiments (i.e. with parameters calibrated for that specific
scenario) to test the proposed metric of human likeness.
Human like behavior can represent a standard for humanoid posture control, and
with the data-set presented in this work, we propose a way to quantify it. Of
course, this does not exclude the possibility that a robot can benefit from
super-human capabilities or peculiar behaviors that are optimized for its non-
human hardware, as motors are obviously different from muscles and almost all
humanoid robots have limb kinematics that is simpler than the corresponding
human one. Comparing FRFs is not the only method available in the literature.
Models can be tested versus human data with statistical tests such as analysis
of variance, focusing on some specific kinematic variables, as presented, for
example, in [Bayon et al., 2020] where the range of motions of the body links
are considered. The direct comparison of FRFs has the main advantage of being
based just on experimental data. In any case, mimicking human behavior may be
intrinsically valuable in human-robot interactions for making movements
understandable and physical interaction natural [Villalobos and Zielińska,
2017]. The presented comparison principle can be extended using several FRFs
obtained with different modalities, e.g. support surface translation and tilt,
and different amplitudes. Performing different tests with the same control
system parameters can highlight some specificity of human behavior such as
nonlinear responses and sensory re-weighting (the importance of specific
sensory inputs can change depending on the conditions). From the practical
point of view, humans can get an advantage from a nonlinear response in the
sense that smaller disturbances can be compensated in a more “relaxed” way
saving energy and effort of the nervous system. This means that the stimulus-
to-sway gain may be relatively larger for smaller stimuli, and smaller for
larger and more dangerous perturbations.
Human likeness is a “human likelihood” once it is defined in terms of the
distance between the tested sample and a distribution of given human
behaviors. In particular the score in eq. 4 resembles a Mahalanobis distance
$D=\sqrt{\Delta^{T}\Sigma^{-1}\Delta}$, but with the addition of weights.
Assuming a joint normal distribution for the weighted FRFs,
$\mathcal{N}(\mu,\Sigma)$, the Mahalanobis distance defines probability
density, where a smaller distance associated with higher probability
density222Normality of variables describing posturography may be implicitly
assumed when the effects on posture are tested through an analysis of
variance, e.g. [Akçay et al., 2021, Jilk et al., 2014].. It comes naturally to
ask how the proposed metric is distributed over the data-set and how it
relates to the likelihood of a sample. Fig 7 shows the distribution of the
human likeness metric over the data-set of human experiments (cumulative
distribution function and relationship between the metric and the Mahalanobis
distance) together with the results on the tested control systems.

Figure 6: FRF obtained with the IC, DEC, and EM as well as the average FRF of
the human data-set. Figure 7: Mahalanobis distance for the samples in the
data-set (red dots) as a function of the metric from eq. 4, and the bootstrap
estimated cumulative distribution of the metric (blue line) and its variance
(light blue bands). The human-likeness score is smaller than the Mahalanobis
distance because if the coefficients in Fig. 2 D that are $\leq 1$, and in
general there is a spread of the Mahalanobis distances of samples associated
with the same D because the weighting almost removes the variation due to
high-frequency components. The three examples are shown for comparison.

The metric we used in this work opens the possibility to make geometrical
considerations for the performance score: i.e. in terms of “How much two
behaviors are different?” or, “Which is the smallest change in the FRF that
would bring the behavior to a certain score?”. This represents an advantage
when one aims to define solutions with better scores.
Translation and tilt of the body support surface can be used in combination.
Exploiting the structure of the PRTS described in §2, two PRTS signals can be
designed with no overlapping peaks $P(f_{Peaks})$ in the power-spectrum. These
two signals can be used together as profiles for the two stimuli and hence two
independent FRF can be produced in a trial. Preliminary studies suggest that,
although human responses to perturbations are in general non-linear, they may
exhibit a linear superposition effect between these two modalities.
Future work will focus on tests with real humanoid robots as soon as possible.
In this work, simulations were used in lieu of real experiments because of the
limited laboratory access due to the pandemic. A performance indicator
evaluating transient response will be developed using a raised cosine velocity
profile instead of the PRTS (that is oriented to the study of steady-state
behavior). This kind of stimulus represents a “softer” and safer version of a
step function input and it was used in previous experiments with human
subjects [Lippi et al., 2014].

##  5 DATA AND CODE

This performance indicator and the data-set used to define the standard for
the benchmark are available through the EUROBENCH project
http://eurobench2020.eu/contact-us/. The EUROBENCH framework includes several
tests and performance indicators provided with the aim of testing the
performance of robots at any stage of development.

## ACKNOWLEDGMENT

This work is supported by the project

COMTEST, a sub-project of EUROBENCH (European Robotic Framework for Bipedal
Locomotion Benchmarking, www.eurobench2020.eu) funded by H2020 Topic ICT
27-2017 under grant agreement number 779963.

## REFERENCES

  * Abel et al., 2020 Abel, M., Kuz, S., Patel, H. J., Petruck, H., Schlick, C. M., Pellicano, A., and Binkofski, F. C. (2020).  Gender effects in observation of robotic and humanoid actions.  Frontiers in psychology, 11:797.
  * Akçay et al., 2021 Akçay, M. E., Lippi, V., and Mergner, T. (2021).  Visual Modulation of Human Responses to Support Surface Translation.  Frontiers in Human Neuroscience, 15:98.
  * AlBakri, 2008 AlBakri, M. (2008).  Development of a mathematical model and simulation environment for the postural robot (PostuRob II).  Technical report, Uniklinik Freiburg.
  * Alexandrov et al., 2005 Alexandrov, A. V., Frolov, A. A., Horak, F. B., Carlson-Kuhta, P., and Park, S. (2005).  Feedback equilibrium control during human standing.  Biological cybernetics, 93(5):309–322.
  * Alexandrov et al., 2017 Alexandrov, A. V., Lippi, V., Mergner, T., Frolov, A. A., Hettich, G., and Husek, D. (2017).  Human-inspired eigenmovement concept provides coupling-free sensorimotor control in humanoid robot.  Front. Neurorobot., 11:22.
  * Assländer et al., 2015 Assländer, L., Hettich, G., and Mergner, T. (2015).  Visual contribution to human standing balance during support surface tilts.  Hum. Mov. Sci., 41:147–164.
  * Atkeson et al., 2015 Atkeson, C. G., Babu, B. P. W., Banerjee, N., Berenson, D., Bove, C. P., Cui, X., DeDonato, M., Du, R., Feng, S., Franklin, P., et al. (2015).  No falls, no resets: Reliable humanoid behavior in the darpa robotics challenge.  In 2015 IEEE-RAS 15th International Conference on Humanoid Robots (Humanoids), pages 623–630. IEEE.
  * Atkeson et al., 2018 Atkeson, C. G., Benzun, P. B., Banerjee, N., Berenson, D., Bove, C. P., Cui, X., DeDonato, M., Du, R., Feng, S., Franklin, P., et al. (2018).  What happened at the DARPA Robotics Challenge finals.  In The DARPA Robotics Challenge Finals: Humanoid Robots to the Rescue, pages 667–684. Springer.
  * Bayon et al., 2020 Bayon, C., Emmens, A., Afschrift, M., Van Wouwe, T., Keemink, A., Van Der Kooij, H., and Van Asseldonk, E. (2020).  Can momentum-based control predict human balance recovery strategies?  IEEE transactions on neural systems and rehabilitation engineering, 28(9):2015–2024.
  * Colasanto et al., 2015 Colasanto, L., Van der Noot, N., and Ijspeert, A. J. (2015).  Bio-inspired walking for humanoid robots using feet with human-like compliance and neuromuscular control.  In 2015 IEEE-RAS 15th international conference on humanoid robots (humanoids), pages 26–32. IEEE.
  * Conti et al., 2018 Conti, R., Giovacchini, F., Saccares, L., Vitiello, N., Pons, J. L., and Torricelli, D. (2018).  What do people expect from benchmarking of bipedal robots? Preliminary results of the EUROBENCH survey.  In International Symposium on Wearable Robotics, pages 132–136. Springer.
  * Goodworth and Peterka, 2018 Goodworth, A. D. and Peterka, R. J. (2018).  Identifying mechanisms of stance control: a single stimulus multiple output model-fit approach.  Journal of Neuroscience Methods, 296:44–56.
  * Guizzo and Ackerman, 2015 Guizzo, E. and Ackerman, E. (2015).  The hard lessons of DARPA’s robotics challenge [News].  IEEE Spectrum, 52(8):11–13.
  * Hettich et al., 2014 Hettich, G., Assländer, L., Gollhofer, A., and Mergner, T. (2014).  Human hip—ankle coordination emerging from multisensory feedback control.  Hum Mov Sci, 37:123–146.
  * Hettich et al., 2013 Hettich, G., Lippi, V., and Mergner, T. (2013).  Human-like sensor fusion mechanisms in a postural control robot.  In Londral, A. E., Encarnacao, P., and Pons, J. L., editors, Proceedings of the International Congress on Neurotechnology, Electronics and Informatics. Vilamoura, Portugal, pages 152–160.
  * Hettich et al., 2015 Hettich, G., Lippi, V., and Mergner, T. (2015).  Human-like sensor fusion implemented in the posture control of a bipedal robot.  In Neurotechnology, Electronics, and Informatics, pages 29–45. Springer.
  * Jilk et al., 2014 Jilk, D. J., Safavynia, S. A., and Ting, L. H. (2014).  Contribution of vision to postural behaviors during continuous support-surface translations.  Experimental Brain Research, 232(1):169–180.
  * Kashiri et al., 2016 Kashiri, N., Ajoudani, A., Caldwell, D. G., and Tsagarakis, N. G. (2016).  Evaluation of hip kinematics influence on the performance of a quadrupedal robot leg.  In Proceedings of the 13th International Conference on Informatics in Control, Automation and Robotics - Volume 1: ICINCO,, pages 205–212. INSTICC, SciTePress.
  * Kim et al., 2009 Kim, S., Kim, C., You, B., and Oh, S. (2009).  Stable whole-body motion generation for humanoid robots to imitate human motions.  In 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 2518–2524. IEEE.
  * Lippi, 2018 Lippi, V. (2018).  Prediction in the context of a human-inspired posture control model.  Robotics and Autonomous Systems.
  * Lippi et al., 2020a Lippi, V., Assländer, L., Akcay, E., and Mergner, T. (2020a).  Body sway responses to pseudorandom support surface translations of vestibular loss subjects resemble those of vestibular able subjects.  Neuroscience Letters, 736:135271.
  * Lippi et al., 2020b Lippi, V., Brands, K. G., and Seel, T. (2020b).  Real-time implementation and evaluation of magnetometerless tracking system for human and humanoid posture control benchmarking based on inertial sensors.  In Proceedings of the 17th International Conference on Informatics in Control, Automation and Robotics - Volume 1: ICINCO, pages 675–680. INSTICC, SciTePress.
  * Lippi et al., 2014 Lippi, V., Hettich, G., and Mergner, T. (2014).  Modeling postural control of support surface translations.  In IEEE Humanoids, Workshop on cognition, perception and postural control for humanoids Madrid, Spain.
  * Lippi and Mergner, 2017 Lippi, V. and Mergner, T. (2017).  Human-derived disturbance estimation and compensation (DEC) method lends itself to a modular sensorimotor control in a humanoid robot.  Front. Neurorobot., 11:49.
  * Lippi et al., 2020c Lippi, V., Mergner, T., Maurer, C., and Seel, T. (2020c).  Performance Indicators of Humanoid Posture Control and Balance Inspired by Human Experiments.  In 2020 International Symposium on Wearable Robotics and Rehabilitation (WeRob).
  * Lippi et al., 2019 Lippi, V., Mergner, T., Seel, T., and Maurer, C. (2019).  COMTEST project: A complete modular test stand for human and humanoid posture control and balance.  In 2019 IEEE-RAS 19th International Conference on Humanoid Robots (Humanoids) Toronto, Canada. October 15-17.
  * Lippi et al., 2016 Lippi, V., Mergner, T., Szumowski, M., Zurawska, M. S., and Zielińska, T. (2016).  Human-inspired humanoid balancing and posture control in frontal plane.  In ROMANSY 21-Robot Design, Dynamics and Control: Proceedings of the 21st CISM-IFToMM Symposium, June 20-23, Udine, Italy, volume 569, pages 285–292. Springer.
  * Mergner et al., 2003 Mergner, T., Maurer, C., and Peterka, R. J. (2003).  A multisensory posture control model of human upright stance.  Progress in Brain Research, 142:189–201.
  * Ott et al., 2016 Ott, C., Henze, B., Hettich, G., Seyde, T. N., Roa, M. A., Lippi, V., and Mergner, T. (2016).  Good posture, good balance: comparison of bioinspired and model-based approaches for posture control of humanoid robots.  IEEE Robotics & Automation Magazine, 23(1):22–33.
  * Oztop et al., 2005 Oztop, E., Franklin, D. W., Chaminade, T., and Cheng, G. (2005).  Human–humanoid interaction: is a humanoid robot perceived as a human?  International Journal of Humanoid Robotics, 2(04):537–559.
  * Pasma et al., 2018 Pasma, J. H., Assländer, L., van Kordelaar, J., de Kam, D., Mergner, T., and Schouten, A. C. (2018).  Evidence in support of the independent channel model describing the sensorimotor control of human stance using a humanoid robot.  Frontiers in Computational Neuroscience, 12:13.
  * Peterka, 2002 Peterka, R. (2002).  Sensorimotor integration in human postural control.  Journal of neurophysiology, 88(3):1097–1118.
  * Peterka, 2003 Peterka, R. J. (2003).  Simplifying the complexities of maintaining balance.  IEEE Engineering in Medicine and Biology Magazine, 22(2):63–68.
  * Singh, 2016 Singh, N. K. (2016).  Detection of hesitant dynamic postural control.  In 2016 IEEE 12th International Colloquium on Signal Processing & Its Applications (CSPA), pages 36–40. IEEE.
  * Torricelli et al., 2014 Torricelli, D., Mizanoor, R. S., Gonzalez, J., Lippi, V., Hettich, G., Asslaender, L., Weckx, M., Vanderborght, B., Dosen, S., Sartori, M., et al. (2014).  Benchmarking human-like posture and locomotion of humanoid robots: A preliminary scheme.  In Conference on Biomimetic and Biohybrid Systems, pages 320–331. Springer.
  * Torricelli et al., 2020 Torricelli, D., Mizanoor, R. S., Lippi, V., Weckx, M., Mathijssen, G., Vanderborght, B., Mergner, T., Lefeber, D., and Pons, J. L. (2020).  Benchmarking human likeness of bipedal robot locomotion: State of the art and future trends.  In Metrics of Sensory Motor Coordination and Integration in Robots and Animals, pages 147–166. Springer.
  * Villalobos and Zielińska, 2017 Villalobos, J. and Zielińska, T. (2017).  Study of postural adjustments for humanoidal helpmates.  Journal of Automation Mobile Robotics and Intelligent Systems, 11\.
  * von Zitzewitz et al., 2013 von Zitzewitz, J., Boesch, P. M., Wolf, P., and Riener, R. (2013).  Quantifying the human likeness of a humanoid robot.  International Journal of Social Robotics, 5(2):263–276.
  * Zebenay et al., 2015 Zebenay, M., Lippi, V., and Mergener, T. (2015).  Human-like humanoid robot posture control.  In Proceedings of the 12th International Conference on Informatics in Control, Automation and Robotics - Volume 2: ICINCO, pages 304–309. INSTICC, SciTePress.