# Identification of Gaussian Process State-Space Models with Particle
Stochastic Approximation EM

Roger Frigola  Fredrik Lindsten  Thomas B. Schön  Carl E. Rasmussen  Dept. of
Engineering, University of Cambridge, United Kingdom.  Div. of Automatic
Control, Linköping University, Sweden.  Dept. of Information Technology,
Uppsala University, Sweden.

###### Abstract

Gaussian process state-space models (GP-SSMs) are a very flexible family of
models of nonlinear dynamical systems. They comprise a Bayesian nonparametric
representation of the dynamics of the system and additional (hyper-)parameters
governing the properties of this nonparametric representation. The Bayesian
formalism enables systematic reasoning about the uncertainty in the system
dynamics. We present an approach to maximum likelihood identification of the
parameters in GP-SSMs, while retaining the full nonparametric description of
the dynamics. The method is based on a stochastic approximation version of the
EM algorithm that employs recent developments in particle Markov chain Monte
Carlo for efficient identification.

###### keywords:

System identification, Bayesian, Non-parametric identification, Gaussian
processes

##  1 Introduction

Inspired by recent developments in robotics and machine learning, we aim at
constructing models of nonlinear dynamical systems capable of quantifying the
uncertainty in their predictions. To do so, we use the Bayesian system
identification formalism whereby degrees of uncertainty and belief are
represented using probability distributions (Peterka, 1981). Our goal is to
identify models that provide error-bars to any prediction. If the data is
informative and the system identification unambiguous, the model will report
high confidence and error-bars will be narrow. On the other hand, if
predictions are made in operating regimes that were not present in the data
used for system identification, we expect the error-bars to be larger.

Nonlinear state-space models are a very general and widely used class of
dynamical system models. They allow for modeling of systems based on observed
input-output data through the use of a latent (unobserved) variable, the
_state_ $\mathbf{x}_{t}\in\mathcal{X}\triangleq\mathbb{R}^{n_{x}}$. A
discrete-time state-space model (SSM) can be described by

$\displaystyle\mathbf{x}_{t+1}$
$\displaystyle=f(\mathbf{x}_{t},\mathbf{u}_{t})+\mathbf{v}_{t},$ (1a)
$\displaystyle\mathbf{y}_{t}$
$\displaystyle=g(\mathbf{x}_{t},\mathbf{u}_{t})+\mathbf{e}_{t},$ (1b)

where $\mathbf{y}_{t}$ represents the output signal, $\mathbf{u}_{t}$ is the
input signal, and $\mathbf{v}_{t}$ and $\mathbf{e}_{t}$ denote i.i.d. noises.
The state transition dynamics is described by the nonlinear function $f$
whereas $g$ links the output data at a given time to the latent state and
input at that same time. For convenience, in the following we will not
explicitly represent the inputs in our formulation. When available, inputs can
be straightforwardly added as additional arguments to the functions $f$ and
$g$.

A common approach to system identification with nonlinear state-space models
consists in defining a parametric form for the functions $f$ and $g$ and
finding the value of the parameters that minimizes a cost function, e.g. the
negative likelihood. Those parametric functions are typically based on
detailed prior knowledge about the system, such as the equations of motion of
an aircraft, or belong to a class of parameterized generic function
approximators, e.g. artificial neural networks (ANNs). In the following, it
will be assumed that no detailed prior knowledge of the system is available to
create a parametric model that adequately captures the complexity of the
dynamical system. As a consequence, we will turn to generic function
approximators. Parametrized nonlinear functions such as radial basis functions
or other ANNs suffer from both theoretical and practical problems. For
instance, a practitioner needs to select a parametric structure for the model,
such as the number of layers and the number of neurons per layer in a neural
network, which are difficult to choose when little is known about the system
at hand. On a theoretical level, fixing the number of parameters effectively
bounds the complexity of the functions that can be fitted to the data
(Ghahramani, 2012). In order to palliate those problems, we will use Gaussian
processes (Rasmussen and Williams, 2006) which provide a practical framework
for Bayesian nonparametric nonlinear system identification.

Gaussian processes (GPs) can be used to identify nonlinear state-space models
by placing GP priors on the unknown functions. This gives rise to the Gaussian
Process State-Space Model (GP-SSM) (Turner et al., 2010; Frigola et al., 2013)
which will be introduced in Section 2. The GP-SSM is a nonparametric model,
though, the GP is in general governed by a (typically) small number of hyper-
parameters, effectively rendering the model semiparametric. In this work, the
hyper-parameters of the model will be estimated by maximum likelihood, while
retaining the full nonparametric richness of the system model. This is
accomplished by analytically marginalizing out the nonparametric part of the
model and using the particle stochastic approximation EM (PSAEM) algorithm by
Lindsten (2013) for estimating the parameters.

Prior work on GP-SSMs includes (Turner et al., 2010), which presented an
approach to maximum likelihood estimation in GP-SSMs based on analytical
approximations and the parameterization of GPs with a pseudo data set. Ko and
Fox (2011) proposed an algorithm to learn (i.e. identify) GP-SSMs based on
observed data which also used weak labels of the unobserved state trajectory.
Frigola et al. (2013) proposed the use of particle Markov chain Monte Carlo to
provide a fully Bayesian solution to the identification of GP-SSMs that did
not need a pseudo data set or weak labels about unobserved states. However,
the fully Bayesian solution requires priors on the model parameters which are
unnecessary when seeking a maximum likelihood solution. Approaches for
filtering and smoothing using already identified GP-SSMs have also been
developed (Deisenroth et al., 2012; Deisenroth and Mohamed, 2012).

##  2 Gaussian Process State-Space Models

###  2.1 Gaussian Processes

Whenever there is an unknown function, GPs allow us to perform Bayesian
inference directly in the space of functions rather than having to define a
parameterized family of functions and perform inference in its parameter
space. GPs can be used as priors over functions that encode vague assumptions
such as smoothness or stationarity. Those assumptions are often less
restrictive than postulating a parametric family of functions.

Formally, a GP is defined as a collection of random variables, any finite
number of which have a joint Gaussian distribution. A GP
$f(\mathbf{x})\in\mathbb{R}$ can be written as

$f(\mathbf{x})\sim\mathcal{GP}\big{(}m(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime})\big{)},$
(2)

where the mean function $m(\mathbf{x})$ and the covariance function
$k(\mathbf{x},\mathbf{x}^{\prime})$ are defined as

$\displaystyle m(\mathbf{x})$ $\displaystyle=\mathbb{E}[f(\mathbf{x})],$ (3a)
$\displaystyle k(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=\mathbb{E}[(f(\mathbf{x})-m(\mathbf{x}))(f(\mathbf{x}^{\prime})-m(\mathbf{x}^{\prime}))].$
(3b)

A finite number of variables from a Gaussian process follow a jointly Gaussian
distribution

$\begin{bmatrix}f(\mathbf{x}_{1})\\\ f(\mathbf{x}_{2})\\\ \vdots\\\
\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m(\mathbf{x}_{1})\\\
m(\mathbf{x}_{2})\\\ \vdots\\\
\end{bmatrix},\begin{bmatrix}k(\mathbf{x}_{1},\mathbf{x}_{1})&k(\mathbf{x}_{1},\mathbf{x}_{2})&\\\
k(\mathbf{x}_{2},\mathbf{x}_{1})&k(\mathbf{x}_{2},\mathbf{x}_{2})&\\\
&&\ddots\\\ \end{bmatrix}\right).$ (4)

We refer the reader to (Rasmussen and Williams, 2006) for a thorough
exposition of GPs.

###  2.2 Gaussian Process State-Space Models

In this article we will focus on problems where there is very little
information about the nature of the state transition function
$f(\mathbf{x}_{t})$ and a GP is used to model it. However, we will consider
that more information is available about $g(\mathbf{x}_{t})$ and hence it will
be modeled by a parametric function. This is reasonable in many cases where
the mapping from states to observations is known, at least up to some
parameters.

The generative probabilistic model for the GP-SSM is fully specified by

$\displaystyle\mathbf{f}_{t+1}\mid\mathbf{x}_{t}$
$\displaystyle\sim\mathcal{GP}\big{(}m_{\boldsymbol{\theta}_{\mathbf{x}}}(\mathbf{x}_{t}),k_{\boldsymbol{\theta}_{\mathbf{x}}}(\mathbf{x}_{t},\mathbf{x}_{t}^{\prime})\big{)},$
(5a) $\displaystyle\mathbf{x}_{t+1}\mid\mathbf{f}_{t+1}$
$\displaystyle\sim\mathcal{N}(\mathbf{x}_{t+1}\mid\mathbf{f}_{t+1},\mathbf{Q}),$
(5b) $\displaystyle\mathbf{y}_{t}\mid\mathbf{x}_{t}$ $\displaystyle\sim
p(\mathbf{y}_{t}\mid\mathbf{x}_{t},\boldsymbol{\theta}_{\mathbf{y}}),$ (5c)

where $\mathbf{f}_{t+1}=f(\mathbf{x}_{t})$ is the value taken by the state
$\mathbf{x}_{t+1}$ after passing through the transition function, but before
the application of process noise $\mathbf{v}_{t+1}$. The Gaussian process in
(5a) describes the prior distribution over the transition function. The GP is
fully specified by its mean function
$m_{\boldsymbol{\theta}_{\mathbf{x}}}(\mathbf{x})$ and its covariance function
$k_{\boldsymbol{\theta}_{\mathbf{x}}}(\mathbf{x}_{t},\mathbf{x}_{t}^{\prime})$,
which are parameterized by the vector of hyper-parameters
$\boldsymbol{\theta}_{\mathbf{x}}$. Equation (5b) describes the addition of
process noise following a zero-mean Gaussian distribution of covariance
$\mathbf{Q}$. We will place no restrictions on the likelihood distribution
(5c) which will be parameterized by a finite-dimensional vector
$\boldsymbol{\theta}_{\mathbf{y}}$. For notational convenience we group all
the (hyper-)parameters into a single vector
$\boldsymbol{\theta}=\\{\boldsymbol{\theta}_{\mathbf{x}},\boldsymbol{\theta}_{\mathbf{y}},\mathbf{Q}\\}$.

##  3 Maximum Likelihood in the GP-SSM

Maximum likelihood (ML) is a widely used frequentist estimator of the
parameters in a statistical model. The ML estimator
$\widehat{\boldsymbol{\theta}}^{\text{ML}}$ is defined as the value of the
parameters that makes the available observations $\mathbf{y}_{0:T}$ as likely
as possible according ot,

$\widehat{\boldsymbol{\theta}}^{\text{ML}}=\operatorname*{arg\,max}_{\boldsymbol{\theta}}\
p(\mathbf{y}_{0:T}\mid\boldsymbol{\theta}).$ (6)

The GP-SSM has two types of latent variables that need to be marginalized
(integrated out) in order to compute the likelihood

$\displaystyle p(\mathbf{y}_{0:T}\mid\boldsymbol{\theta})=\int
p(\mathbf{y}_{0:T},\mathbf{x}_{0:T},\mathbf{f}_{1:T}\mid\boldsymbol{\theta})\
\text{d}\mathbf{x}_{0:T}\,\text{d}\mathbf{f}_{1:T}$ $\displaystyle=\int
p(\mathbf{y}_{0:T}\mid\mathbf{x}_{0:T},\boldsymbol{\theta})\left(\int
p(\mathbf{x}_{0:T},\mathbf{f}_{1:T}\mid\boldsymbol{\theta})\
\text{d}\mathbf{f}_{1:T}\right)\text{d}\mathbf{x}_{0:T}.$ (7)

Following results from Frigola et al. (2013), the latent variables
$\mathbf{f}_{1:T}$ can be marginalized analytically. This is equivalent to
integrating out the uncertainty in the unknown function $f$ and working
directly with a prior over the state trajectories
$p(\mathbf{x}_{0:T}\mid\boldsymbol{\theta})$ that encodes the assumptions
(e.g. smoothness) of $f$ specified in (5a). The prior over trajectories can be
factorized as

$\displaystyle
p(\mathbf{x}_{0:T}\mid\boldsymbol{\theta})=p(\mathbf{x}_{0}\mid\boldsymbol{\theta})\prod_{t=1}^{T}p(\mathbf{x}_{t}\mid\boldsymbol{\theta},\mathbf{x}_{0:t-1}).$
(8)

Using standard expressions for GP prediction, the one-step predictive density
is given by

$\displaystyle
p(\mathbf{x}_{t}\mid\boldsymbol{\theta},\mathbf{x}_{0:t-1})=\mathcal{N}\big{(}\mathbf{x}_{t}\mid\boldsymbol{\mu}_{t}(\mathbf{x}_{0:t-1}),\mathbf{\Sigma}_{t}(\mathbf{x}_{0:t-1})\big{)},$
(9a) where $\displaystyle\boldsymbol{\mu}_{t}(\mathbf{x}_{0:t-1})$
$\displaystyle=\mathbf{m}_{t-1}+\mathbf{K}_{t-1,0:t-2}\widetilde{\mathbf{K}}^{-1}_{0:t-2}\
(\mathbf{x}_{1:t-1}-\mathbf{m}_{0:t-2}),$ (9b)
$\displaystyle\mathbf{\Sigma}_{t}(\mathbf{x}_{0:t-1})$
$\displaystyle=\widetilde{\mathbf{K}}_{t-1}-\mathbf{K}_{t-1,0:t-2}\widetilde{\mathbf{K}}^{-1}_{0:t-2}\mathbf{K}_{t-1,0:t-2}^{\top},$
(9c)

for $t\geq 2$ and $\boldsymbol{\mu}_{1}(\mathbf{x}_{0})=\mathbf{m}_{0}$,
$\mathbf{\Sigma}_{1}(\mathbf{x}_{0})=\widetilde{\mathbf{K}}_{0}$. Here we have
defined the mean vector
$\mathbf{m}_{0:t-1}\triangleq\begin{bmatrix}m({\mathbf{x}}_{0})^{\top}&\dots&m({\mathbf{x}}_{t-1})^{\top}\end{bmatrix}^{\top}$
and the $(n_{x}t)\times(n_{x}t)$ positive definite matrix $\mathbf{K}_{0:t-1}$
with block entries
$[\mathbf{K}_{0:t-1}]_{i,j}=k({\mathbf{x}}_{i-1},{\mathbf{x}}_{j-1})$. These
matrices use two sets of indices, as in $\mathbf{K}_{t-1,0:t-2}$, to refer to
the off-diagonal blocks of $\mathbf{K}_{0:t-1}$. We also define
$\widetilde{\mathbf{K}}_{0:t-1}=\mathbf{K}_{0:t-1}+\mathbf{I}_{t}\otimes\mathbf{Q}$,
where $\otimes$ denotes the Kronecker product.

Using (8) we can thus write the likelihood (7) as

$\displaystyle p(\mathbf{y}_{0:T}\mid\boldsymbol{\theta})=\int
p(\mathbf{y}_{0:T}\mid\mathbf{x}_{0:T},\boldsymbol{\theta})p(\mathbf{x}_{0:T}\mid\boldsymbol{\theta})\
\text{d}\mathbf{x}_{0:T}.$ (10)

The integration with respect to $\mathbf{x}_{0:T}$, however, is not
analytically tractable. This difficulty will be addressed in the subsequent
section.

A GP-SSM can be seen as a hierarchical probabilistic model which describes a
prior over the latent state trajectories
$p(\mathbf{x}_{0:T}\mid\boldsymbol{\theta}_{\mathbf{x}},\mathbf{Q})$ and links
this prior with the observed data via the likelihood
$p(\mathbf{y}_{t}\mid\mathbf{x}_{t},\boldsymbol{\theta}_{\mathbf{y}})$. Direct
application of maximum likelihood on
$p(\mathbf{y}_{t}\mid\mathbf{x}_{t},\boldsymbol{\theta}_{\mathbf{y}})$ to
obtain estimates of the state trajectory and likelihood parameters would
invariably result in over-fitting. However, by introducing a prior on the
state trajectories111A prior over the state trajectories is not an exclusive
feature of GP-SSMs. Linear-Gaussian state-space models, for instance, also
describe a prior distribution over state trajectories:
$p(\mathbf{x}_{1:T}\mid\mathbf{A},\mathbf{B},\mathbf{Q},\mathbf{x}_{0})$. and
marginalizing them as in (10), we obtain the so-called marginal likelihood.
Maximization of the marginal likelihood with respect to the parameters results
in a procedure known as type II maximum likelihood or empirical Bayes (Bishop,
2006). Empirical Bayes reduces the risk of over-fitting since it automatically
incorporates a trade-off between model fit and model complexity, a property
often known as Bayesian Occam’s razor (Ghahramani, 2012).

##  4 Particle Stochastic Approximation EM

As pointed out above, direct evaluation of the likelihood (10) is not possible
for a GP-SSM. However, by viewing the latent states $\mathbf{x}_{0:T}$ as
missing data, we are able to evaluate the _complete data_ log-likelihood

$\log p(\mathbf{y}_{0:T},\mathbf{x}_{0:T}\mid\boldsymbol{\theta})=\log
p(\mathbf{y}_{0:T}\mid\mathbf{x}_{0:T},\boldsymbol{\theta})+\log
p(\mathbf{x}_{0:T}\mid\boldsymbol{\theta}),$ (11)

by using (5c) and (9). We therefore turn to the Expectation Maximization (EM)
algorithm (Dempster et al., 1977). The EM algorithm uses (11) to construct a
surrogate cost function for the ML problem, defined as

$\displaystyle Q({}$
$\displaystyle\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})=\mathbb{E}_{\boldsymbol{\theta}^{\prime}}[\log
p(\mathbf{y}_{0:T},\mathbf{x}_{0:T}\mid\boldsymbol{\theta})\mid y_{0:T}]$
$\displaystyle=\int\log
p(\mathbf{y}_{0:T},\mathbf{x}_{0:T}\mid\boldsymbol{\theta})p(\mathbf{x}_{0:T}\mid\mathbf{y}_{0:T},\boldsymbol{\theta}^{\prime})\text{d}\mathbf{x}_{0:T}.$
(12)

It is an iterative procedure that maximizes (10) by iterating two steps,
expectation (E) and maximization (M),

  * (E)

Compute $Q(\boldsymbol{\theta},\boldsymbol{\theta}_{k-1})$.

  * (M)

Compute
$\boldsymbol{\theta}_{k}=\operatorname*{arg\,max}_{\boldsymbol{\theta}}Q(\boldsymbol{\theta},\boldsymbol{\theta}_{k-1})$.

The resulting sequence $\\{\boldsymbol{\theta}_{k}\\}_{k\geq 0}$ will, under
weak assumptions, converge to a stationary point of the likelihood
$p(\mathbf{y}_{0:T}\mid\boldsymbol{\theta})$.

To implement the above procedure we need to compute the integral in (12),
which in general is not computationally tractable for a GP-SSM. To deal with
this difficulty, we employ a Monte-Carlo-based implementation of the EM
algorithm, referred to as PSAEM (Lindsten, 2013). This procedure is a
combination of stochastic approximation EM (SAEM) (Delyon et al., 1999) and
particle Markov chain Monte Carlo (PMCMC) (Andrieu et al., 2010; Lindsten et
al., 2012). As illustrated by Lindsten (2013), PSAEM is a competitive
alternative to particle-smoothing-based EM algorithms (e.g. (Schön et al.,
2011; Olsson et al., 2008)), as it enjoys better convergence properties and
has a much lower computational cost. The method maintains a stochastic
approximation of the auxiliary quantity (12),
$\widehat{Q}_{k}(\boldsymbol{\theta})\approx
Q(\boldsymbol{\theta},\boldsymbol{\theta}_{k-1})$. This approximation is
updated according to

$\displaystyle\widehat{Q}_{k}(\boldsymbol{\theta})=(1-\gamma_{k})\widehat{Q}_{k-1}(\boldsymbol{\theta})+\gamma_{k}\log
p(\mathbf{y}_{0:T},\mathbf{x}_{0:T}[k]\mid\boldsymbol{\theta}).$ (13)

Here, $\\{\gamma_{k}\\}_{k\geq 0}$ is a sequence of step sizes, satisfying the
usual stochastic approximation conditions: $\sum_{k}\gamma_{k}=\infty$ and
$\sum_{k}\gamma_{k}^{2}<\infty$. A typical choice is to take
$\gamma_{k}=k^{-p}$ with $p\in\,]0.5,1]$, where a smaller value of $p$ gives a
more rapid convergence at the cost of higher variance. In the vanilla SAEM
algorithm, $\mathbf{x}_{0:T}[k]$ is a draw from the smoothing distribution
$p(\mathbf{x}_{0:T}\mid\mathbf{y}_{0:T},\boldsymbol{\theta}_{k-1})$. In this
setting, Delyon et al. (1999) show that using the stochastic approximation
(13) instead of (12) in the EM algorithm results in a valid method, i.e.
$\\{\boldsymbol{\theta}_{k}\\}_{k\geq 0}$ will still converge to a maximizer
of $p(\mathbf{y}_{0:T}\mid\boldsymbol{\theta})$.

The PSAEM algorithm is an extension of SAEM, which is useful when it is not
possible to sample directly from the joint smoothing distribution. This is
indeed the case in our setting. Instead of sampling from the smoothing
distribution, the sample trajectory $\mathbf{x}_{0:T}[k]$ in (13) may be drawn
from an ergodic Markov kernel, leaving the smoothing distribution invariant.
Under suitable conditions on the kernel, this will not violate the validity of
SAEM; see (Andrieu and Vihola, 2011; Andrieu et al., 2005).

In PSAEM, this Markov kernel on the space of trajectories, denoted as
${P_{\boldsymbol{\theta}}^{N}(\mathbf{x}_{0:T}^{\star}\mid\mathbf{x}_{0:T}^{\prime})}$,
is constructed using PMCMC theory. In particular, we use the method by
Lindsten et al. (2012), particle Gibbs with ancestor sampling (PGAS). We have
previously used PGAS for Bayesian identification of GP-SSMs (Frigola et al.,
2013). PGAS is a sequential Monte Carlo method, akin to a standard particle
filter (see e.g. (Doucet and Johansen, 2011; Gustafsson, 2010)), but with the
difference that one particle at each time point is specified _a priori_. These
reference states, denoted as $\mathbf{x}_{0:T}^{\prime}$, can be thought of as
guiding the particles of the particle filter to the “correct” regions of the
state-space. More formally, as shown by Lindsten et al. (2012), PGAS defines a
Markov kernel which leaves the joint smoothing distribution invariant, i.e.
for any $\boldsymbol{\theta}$,

$\displaystyle\int
P_{\boldsymbol{\theta}}^{N}(\mathbf{x}_{0:T}^{\star}\mid\mathbf{x}_{0:T}^{\prime})p(\mathbf{x}_{0:T}^{\prime}\mid\mathbf{y}_{0:T},\boldsymbol{\theta})\text{d}\mathbf{x}_{0:T}^{\prime}=p(\mathbf{x}_{0:T}^{\star}\mid\mathbf{y}_{0:T},\boldsymbol{\theta}).$
(14)

The PGAS kernel is indexed by $N$, which is the number of particles used in
the underlying particle filter. Note in particular that the desired property
(14) holds for any ${N\geq 1}$, i.e. the number of particles only affects the
mixing of the Markov kernel. A larger $N$ implies faster mixing, which in turn
results in better approximations of the auxiliary quantity (13). However, it
has been experienced in practice that the correlation between consecutive
trajectories drops of quickly as $N$ increases (Lindsten et al., 2012;
Lindsten and Schön, 2013), and for many models a moderate $N$ (e.g. in the
range 5–20) is enough to get a rapidly mixing kernel. We refer to (Lindsten,
2013; Lindsten et al., 2012) for details. We conclude by noting that it is
possible to generate a sample $\mathbf{x}_{0:T}[k]\sim
P_{\boldsymbol{\theta}[k-1]}^{N}(\,\cdot\mid\mathbf{x}_{0:T}[k-1])$ by running
a particle-filter-like algorithm. This method is given as Algorithm 1 in
(Lindsten et al., 2012) and is described specifically for GP-SSMs in Section 3
of (Frigola et al., 2013).

Next, we address the M-step of the EM algorithm. Maximizing the quantity (13)
will typically not be possible in closed form. Instead, we make use of a
numerical optimization routine implementing a quasi-Newton method (BFGS).
Using (11), the gradient of the complete data log-likelihood can be written as

$\displaystyle\frac{\partial}{\partial\boldsymbol{\theta}}\log
p(\mathbf{y}_{0:T},\mathbf{x}_{0:T}\mid\boldsymbol{\theta})=\sum_{t=0}^{T}\frac{\partial}{\partial\boldsymbol{\theta}}\log
p(\mathbf{y}_{t}\mid\mathbf{x}_{t},\boldsymbol{\theta})$ $\displaystyle\hskip
10.00002pt+\sum_{t=1}^{T}\frac{\partial}{\partial\boldsymbol{\theta}}\log
p(\mathbf{x}_{t}\mid\mathbf{x}_{0:t-1},\boldsymbol{\theta})+\frac{\partial}{\partial\boldsymbol{\theta}}\log
p(\mathbf{x}_{0}\mid\boldsymbol{\theta}),$ (15)

where the individual terms can be computed using (5c) and (9), respectively.
The resulting PSAEM algorithm for learning of GP-SSMs is summarized in
Algorithm 1.

Algorithm 1 PSAEM for GP-SSMs

  1. 1.

Set $\boldsymbol{\theta}_{0}$ and $\mathbf{x}_{0:T}[0]$ arbitrarily. Set
$\widehat{Q}_{0}(\boldsymbol{\theta})\equiv 0$.

  2. 2.

For $k\geq 1$:

    1. (a)

Simulate $\mathbf{x}_{0:T}[k]\sim
P_{\boldsymbol{\theta}[k-1]}^{N}(\,\cdot\mid\mathbf{x}_{0:T}[k-1])$ (run
Algorithm 1 in (Lindsten et al., 2012) and set $\mathbf{x}_{0:T}[k]$ to one of
the particle trajectories with probabilities given by their importance
weights).

    2. (b)

Update $\widehat{Q}_{k}(\boldsymbol{\theta})$ according to (13).

    3. (c)

Compute
$\boldsymbol{\theta}_{k}=\operatorname*{arg\,max}_{\boldsymbol{\theta}}\widehat{Q}_{k}(\boldsymbol{\theta})$.

A particular feature of the proposed approach is that it performs smoothing
even when the state transition function is not yet explicitly defined. Once
samples from the smoothing distribution have been obtained it is then possible
to analytically describe the state transition probability density (see
(Frigola et al., 2013) for details). This contrasts with the standard
procedure where the smoothing distribution is found using a given state
transition density.

##  5 Experimental Results

In this section we present the results of applying PSAEM to identify various
dynamical systems.

###  5.1 Identification of a Linear System

Although GP-SSMs are particularly suited to nonlinear system identification,
we start by illustrating their behavior when identifying the following linear
system

$\displaystyle\mathbf{x}_{t+1}$
$\displaystyle=0.8\,\mathbf{x}_{t}+3\,\mathbf{u}_{t}+\mathbf{v}_{t},$
$\displaystyle\mathbf{v}_{t}\sim\mathcal{N}(0,1.5),$ (16a)
$\displaystyle\mathbf{y}_{t}$
$\displaystyle=2\,\mathbf{x}_{t}+\mathbf{e}_{t},$
$\displaystyle\mathbf{e}_{t}\sim\mathcal{N}(0,1.5),$ (16b)

excited by a periodic input. The GP-SSM can model this linear system by using
a linear covariance function for the GP. This covariance function imposes, in
a somehow indirect fashion, that the state-transition function in (16a) must
be linear. A GP-SSM with linear covariance function is formally equivalent to
a linear state-space model where a Gaussian prior is placed over the, unknown
to us, parameters ($A=0.8$ and $B=3$) (Rasmussen and Williams, 2006, Section
2.1). The hyper-parameters of the covariance function are equivalent to the
variances of a zero-mean prior over $A$ and $B$. Therefore, the application of
PSAEM to this particular GP-SSM can be interpreted as finding the hyper-
parameters of a Gaussian prior over the parameters of the linear model that
maximize the likelihood of the observed data whilst marginalizing over $A$ and
$B$. In addition, the likelihood will be simultaneously optimized with respect
to the process noise and measurement noise variances ($q$ and $r$
respectively).

Figure 1: Convergence of parameters when learning a linear system using a
linear covariance function. Figure 2: Linear dynamical system learned using a
GP-SSM with linear covariance function. Predictions (a) on training data, and
(b) on test data (see text for more details).

Figure 1 shows the convergence of the GP hyper-parameters ($l_{x}$ and
$l_{u}$) and noise parameters with respect to the PSAEM iteration. In order to
judge the quality of the learned GP-SSM we evaluate its predictive performance
on the data set used for learning (training set) and on an independent data
set generated from the same dynamical system (test set). The GP-SSM can make
probabilistic predictions which report the uncertainty arising from the fact
that only a finite amount of data is observed.

Figure 2 displays the predicted value of $\mathbf{f}_{t+1}-\mathbf{x}_{t}$
versus the true value. Recall that $\mathbf{f}_{t+1}-\mathbf{x}_{t}$ is
equivalent to the step taken by the state in one single transition before
process noise is added: $f(\mathbf{x}_{t},\mathbf{u}_{t})-\mathbf{x}_{t}$. One
standard deviation error bars from the predictive distribution have also been
plotted. Perfect predictions would lie on the unit slope line. We note that
although the predictions are not perfect, error-bars tend to be large in
predictions that are far from the true value and narrower for predictions that
are closer to the truth. This is the desired outcome since the goal of the GP-
SSM is to represent the uncertainty in its predictions.

Figure 3: Convergence of parameters when learning a linear system using a
squared exponential covariance function. Figure 4: Linear dynamical system
learned using a GP-SSM with squared exponential covariance function.
Predictions (a) on training data, and (b) on test data.

We now move into a scenario in which the data is still generated by the linear
dynamical system in (16) but we pretend that we are not aware of its
linearity. In this case, a covariance function able to model nonlinear
transition functions is a judicious choice. We use the squared exponential
covariance function which imposes the assumption that the state transition
function is smooth and infinitely differentiable (Rasmussen and Williams,
2006). Figure 3 shows, for a PSAEM run, the convergence of the covariance
function hyper-parameters (length-scales $\lambda_{x}$ and $\lambda_{u}$ and
signal variance $\sigma_{f}$) and also the convergence of the noise
parameters.

The predictive performance on training data and independent test data is
presented in Figure 4. Interestingly, in the panel corresponding to training
data (a), there is particularly poor prediction that largely underestimates
the value of the state transition. However, the variance for this prediction
is very high which indicates that the identified model has little confidence
in it. In this particular case, the mean of the prediction is 2.5 standard
deviations away from the true value of the state transition.

###  5.2 Identification of a Nonlinear System

Figure 5: Nonlinear dynamical system with one state and one input. The black
mesh represents the ground truth dynamics function and the colored surface is
the mean of the identified function. Color is proportional to the standard
deviation of the identified function (red represents high uncertainty and blue
low uncertainty). Figure 6: State trajectory from a test data set (solid black
line). One step ahead predictions made with the identified model are depicted
by a dashed line (mean) and a colored interval at $\pm 1$ standard deviation
(including process noise).

GP-SSMs are particularly powerful for nonlinear system identification when it
is not possible to create a parametric model of the system based on detailed
knowledge about its dynamics. To illustrate this capability of GP-SSMs we
consider the nonlinear dynamical system

$\displaystyle\mathbf{x}_{t+1}$
$\displaystyle=a\mathbf{x}_{t}+b\frac{\mathbf{x}_{t}}{1+\mathbf{x}_{t}^{2}}+c\mathbf{u}_{t}+\mathbf{v}_{t},$
$\displaystyle\mathbf{v}_{t}\sim\mathcal{N}(0,q),$ (17a)
$\displaystyle\mathbf{y}_{t}$
$\displaystyle=d\mathbf{x}_{t}^{2}+\mathbf{e}_{t},$
$\displaystyle\mathbf{e}_{t}\sim\mathcal{N}(0,r),$ (17b)

with parameters $(a,b,c,d,q,r)=(0.5,25,8,0.05,10,1)$ and a known input
$\mathbf{u}_{t}=\cos(1.2(t+1))$. One of the challenging properties of this
system is that the quadratic measurement function (17b) tends to induce a
bimodal distribution in the marginal smoothing distribution. For instance, if
we were to consider only one measurement in isolation and $r=0$ we would have
$\mathbf{x}_{t}=\pm\sqrt{\frac{\mathbf{y}_{t}}{d}}$. Moreover, the state
transition function (17a) exhibits a very sharp gradient in the
$\mathbf{x}_{t}$ direction at the origin, but is otherwise parsimonious as
$\mathbf{x}_{t}\rightarrow\pm\infty$.

Again, we pretend that detailed knowledge about the particular form of (17a)
is not available to us. We select a covariance function that consists of a
Matérn covariance function in the $\mathbf{x}$ direction and a squared
exponential in the $\mathbf{u}$ direction. The Matérn covariance function
imposes less smoothness constraints than the squared exponential (Rasmussen
and Williams, 2006) and is therefore more suited to model functions that can
have sharp transitions.

Figure 5 shows the true state transition dynamics function (black mesh) and
the identified function as a colored surface. Since the identified function
from the GP-SSM comes in the form of a probability distribution over
functions, the surface is plotted at
$\mathbb{E}[\mathbf{f}^{*}|\mathbf{x}^{*},\mathbf{u}^{*},\mathbf{y}_{0:T}]$
where the symbol ∗ denotes test points. The standard deviation of
$\mathbf{f}^{*}$, which represents our uncertainty about the actual value of
the function, is depicted by the color of the surface. Figure 6 shows the one
step ahead predictive distributions
$p(\mathbf{x}_{t+1}^{*}|\mathbf{x}_{t}^{*},\mathbf{u}_{t}^{*},\mathbf{y}_{0:T})$
on a test data set.

##  6 Conclusions

GP-SSMs allow for a high degree of flexibility when addressing the nonlinear
system identification problem by making use of Bayesian nonparametric system
models. These models enable the incorporation of high-level assumptions, such
as smoothness of the transition function, while still being able to capture a
wide range of nonlinear dynamical functions. Furthermore, the GP-SSM is
capable of making probabilistic predictions that can be useful in adaptive
control and robotics, when the control strategy might depend on the
uncertainty in the dynamics. Our particle-filter-based maximum likelihood
inference of the model hyper-parameters preserves the full nonparametric
richness of the model. In addition, marginalization of the dynamical function
effectively averages over all possible dynamics consistent with the GP prior
and the data, and hence provides a strong safeguard against overfitting.

## References

  * Andrieu and Vihola (2011) C. Andrieu and M. Vihola.  Markovian stochastic approximation with expanding projections.  arXiv.org, arXiv:1111.5421, November 2011.
  * Andrieu et al. (2005) C. Andrieu, E. Moulines, and P. Priouret.  Stability of stochastic approximation under verifiable conditions.  _SIAM Journal on Control and Optimization_ , 44(1):283–312, 2005.
  * Andrieu et al. (2010) C. Andrieu, A. Doucet, and R. Holenstein.  Particle Markov chain Monte Carlo methods.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 72(3):269–342, 2010.  ISSN 1467-9868.
  * Bishop (2006) C. M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Deisenroth and Mohamed (2012) M. P. Deisenroth and S. Mohamed.  Expectation Propagation in Gaussian process dynamical systems.  In P. Bartlett, F.C.N. Pereira, C.J.C. Burges, L. Bottou, and K.Q. Weinberger, editors, _Advances in Neural Information Processing Systems (NIPS) 25_ , pages 2618–2626. 2012.
  * Deisenroth et al. (2012) M. P. Deisenroth, R. D. Turner, M. F. Huber, U. D. Hanebeck, and C. E. Rasmussen.  Robust filtering and smoothing with Gaussian processes.  _IEEE Transactions on Automatic Control_ , 57(7):1865 –1871, july 2012.  ISSN 0018-9286.  10.1109/TAC.2011.2179426.
  * Delyon et al. (1999) B. Delyon, M. Lavielle, and E. Moulines.  Convergence of a stochastic approximation version of the EM algorithm.  _The Annals of Statistics_ , 27(1):pp. 94–128, 1999.
  * Dempster et al. (1977) A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 39(1):pp. 1–38, 1977.
  * Doucet and Johansen (2011) A. Doucet and A. Johansen.  A tutorial on particle filtering and smoothing: Fifteen years later.  In D. Crisan and B. Rozovskii, editors, _The Oxford Handbook of Nonlinear Filtering_. Oxford University Press, 2011.
  * Frigola et al. (2013) R. Frigola, F. Lindsten, T. B. Schön, and C. E. Rasmussen.  Bayesian inference and learning in Gaussian process state-space models with particle MCMC.  In L. Bottou, C.J.C. Burges, Z. Ghahramani, M. Welling, and K.Q. Weinberger, editors, _Advances in Neural Information Processing Systems (NIPS) 26_. 2013.
  * Ghahramani (2012) Z. Ghahramani.  Bayesian nonparametrics and the probabilistic approach to modelling.  _Philosophical Transactions of the Royal Society A_ , 2012.
  * Gustafsson (2010) F. Gustafsson.  Particle filter theory and practice with positioning applications.  _IEEE Aerospace and Electronic Systems Magazine_ , 25(7):53–82, 2010.
  * Ko and Fox (2011) J. Ko and D. Fox.  Learning GP-BayesFilters via Gaussian process latent variable models.  _Autonomous Robots_ , 30(1):3–23, 2011.
  * Lindsten (2013) F. Lindsten.  An efficient stochastic approximation EM algorithm using conditional particle filters.  In _Proceedings of the 38th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)_ , pages 6274–6278, 2013.
  * Lindsten and Schön (2013) F. Lindsten and T. B. Schön.  Backward simulation methods for Monte Carlo statistical inference.  _Foundations and Trends in Machine Learning_ , 6(1):1–143, 2013.
  * Lindsten et al. (2012) F. Lindsten, M. Jordan, and T. B. Schön.  Ancestor sampling for particle Gibbs.  In P. Bartlett, F.C.N. Pereira, C.J.C. Burges, L. Bottou, and K.Q. Weinberger, editors, _Advances in Neural Information Processing Systems (NIPS) 25_ , pages 2600–2608. 2012.
  * Olsson et al. (2008) J. Olsson, R. Douc, O. Cappé, and E. Moulines.  Sequential Monte Carlo smoothing with application to parameter estimation in nonlinear state-space models.  _Bernoulli_ , 14(1):155–179, 2008.
  * Peterka (1981) V. Peterka.  Bayesian system identification.  _Automatica_ , 17(1):41 – 53, 1981.
  * Rasmussen and Williams (2006) C. E. Rasmussen and C. K. I. Williams.  _Gaussian Processes for Machine Learning_.  MIT Press, 2006.
  * Schön et al. (2011) T. B. Schön, A. Wills, and B. Ninness.  System identification of nonlinear state-space models.  _Automatica_ , 47(1):39–49, 2011.
  * Turner et al. (2010) R. Turner, M. P. Deisenroth, and C. E. Rasmussen.  State-space inference and learning with Gaussian processes.  In Yee Whye Teh and Mike Titterington, editors, _13th International Conference on Artificial Intelligence and Statistics_ , volume 9 of _W &CP_, pages 868–875, Chia Laguna, Sardinia, Italy, May 13–15 2010\.


# Infrared (in)sensitivity of relativistic effects in cosmological observable
statistics

Ermis Mitsou  Center for Theoretical Astrophysics and Cosmology, Institute for
Computational Science, University of Zurich, CH–8057 Zurich, Switzerland
Jaiyul Yoo<EMAIL_ADDRESS>Center for Theoretical Astrophysics and
Cosmology, Institute for Computational Science, University of Zurich, CH–8057
Zurich, Switzerland  Physics Institute, University of Zurich,
Winterthurerstrasse 190, CH–8057, Zurich, Switzerland  Matteo Magi  Center for
Theoretical Astrophysics and Cosmology, Institute for Computational Science,
University of Zurich, CH–8057 Zurich, Switzerland

###### Abstract

The relativistic effects in cosmological observables contain critical
information about the initial conditions and gravity on large scales. Compared
to the matter density fluctuation, some of these relativistic contributions
scale with negative powers of comoving wave number, implying a growing
sensitivity to infrared modes. However, this can be inconsistent with the
equivalence principle and can also lead to infrared divergences in the
observable $N$-point statistics. Recent perturbative calculations have shown
that this infrared sensitivity is indeed spurious due to subtle cancellations
in the cosmological observables that have been missed in the bulk of the
literature. Here we demonstrate that the cosmological observable statistics
are infrared-insensitive in a general and fully non-linear way, assuming
diffeomorphism invariance and adiabatic fluctuations on large scales.

Introduction.— In cosmology, constraining model parameters with observations
requires relating the $N$-point correlation functions of cosmological
observables to the $N$-point correlation functions of space-time fields at
some initial hypersurface (see, e.g., [1, 2, 3, 4]). Cosmological observables
$O$, such as the luminosity distance, galaxy number density or cosmic
microwave background temperature, are typically a function of the observed
redshift $z$ of a source and its direction $\hat{n}$ in the sky, together
parametrizing the past light-cone. The initial condition fields, such as the
comoving-gauge curvature perturbation $\mathcal{R}$ and tensor perturbation
$\gamma_{ij}$ for minimal models, are typically considered in the Fourier
space of the initial hypersurface, parametrized by $k^{i}$, so that
statistical homogeneity and isotropy simplify their statistics.

As an example, the two-point correlation function of the cosmological
observables in linear perturbation theory for the scalar sector reads

$\langle
O(z,\hat{n})\,O(z^{\prime},\hat{n}^{\prime})\rangle=\bar{O}(z)\,\bar{O}(z^{\prime})\left[1+\sum_{l=0}^{\infty}{\cal
L}_{l}(\mu)\,C_{l}(z,z^{\prime})\right]\,,$ (1)

where ${\cal L}_{l}$ are the Legendre polynomials with cosine angle
$\mu:=\hat{n}\cdot\hat{n}^{\prime}$,

$C_{l}(z,z^{\prime})=\int d\ln
k~{}\mathcal{T}_{l}(z;k)\,\mathcal{T}_{l}(z^{\prime};k)\,\Delta^{2}_{\mathcal{R}}(k)~{},$
(2)

is the angular power spectrum of the cosmological observables, and
$\Delta^{2}_{\mathcal{R}}:=k^{3}P_{\mathcal{R}}(k)/2\pi^{2}\propto
k^{n_{s}-1}$ is the dimensionless power spectrum of the primordial
fluctuations $\mathcal{R}$ with spectral index $n_{s}$. All the information
about the evolution of the universe and the light propagation to the observer
is then contained in the kernel or the transfer function
$\mathcal{T}_{l}(z;k)$. The generalization of this relation in Eq. (2) to
higher-order perturbation theory and other $N$-point statistics retains the
same qualitative form of a convolution with kernels that depend on multipole
indices, redshifts and initial-hypersurface Fourier wave-vectors 111The
Fourier coordinates $k^{i}$ considered here should not be confused with the
ones observers use to map the sky, as the two types are equivalent only in the
Newtonian limit. Indeed, the “observer’s Fourier space” is the dual of
Cartesian coordinates inferred from the data $(z,\hat{n})$ through some
fiducial cosmology and thus parametrizing the past light-cone, not a constant-
time hypersurface. Consequently, in this observational Fourier space large
scales correspond to early times, meaning that there is either no well-defined
infrared limit (finite particle horizon), or that limit is trivial since
structure is absent in the distant past..

The Fourier dependence of these kernels can then be organized in decreasing
orders of spatial derivatives, i.e. the derivatives that act on the primordial
fields $\mathcal{R}$ before taking the ensemble average. This ordering makes
particular sense, since large-scale surveys have historically progressed by
probing from small to larger scales (see, e.g., [6, 7]), so higher-order
derivative terms have been most relevant in practical situations. At the level
of the kernels, such as $\mathcal{T}_{l}$ in Eq. (2), this corresponds to
ordering the terms by $({\cal H}/k)^{n\geq 0}$ factors with respect to the
dominating contributions on small scales, or the matter density fluctuation
$\delta_{m}\sim k^{2}\mathcal{R}$, where ${\cal H}$ is the conformal Hubble
parameter. While the lowest-order approximation ${\cal H}/k\to 0$ has been
accurate so far, the upcoming surveys [8, 9, 10, 11, 12, 13] will be sensitive
to the contributions $\sim({\cal H}/k)^{n>0}$. These terms contain information
about the initial conditions in the early Universe and relativistic effects,
both in field equations and the light propagation between source and observer.

If the integrand leading to the observable $N$-point statistics contains low
enough powers of $k$ (or $n>0$), then the result can be dominated by infrared
mode contributions, i.e. with their influence growing indefinitely as $k\to
0$. In general the observer-source system is in free-fall and of finite
extent, so it should gradually lose sensitivity to fluctuations of wavelength
longer than the separation between the source and the observer, because the
infrared limit that is a uniform gravitational potential and force is
equivalent to an accelerated coordinate system (equivalence principle). In
fact, in the most extreme case the Fourier integral in Eq. (2) can even be
infrared-divergent, which is even more unphysical for observations involving a
finite space-time patch. As it turns out [14, 15, 16, 17, 18, 19, 20, 21, 22,
23, 24, 25, 26, 27, 28, 29], the cosmological observables of most interest do
contain relativistic effects of low enough $k$-powers ($n>0$) for this to
happen. Importantly, the current standard approach to modeling the
cosmological observables in the literature leads to infrared-divergences. As a
result, the bulk of the works in the literature offers an inaccurate modeling
of the infrared contributions to the observable $N$-point statistics on large
scales, precisely the region of interest for upcoming observations.

In this Letter we provide a general and fully non-linear proof that the
statistics of cosmological observables remain insensitive to infrared modes
when relativistic effects are taken into account, assuming adiabatic
fluctuations on large enough scales. If one expands the full $k$-dependence of
the kernels, such as $\mathcal{T}_{l}$ in Eq. (2), then the lowest powers of
$k$ from individual relativistic effects cancel each other out exactly. The
precise total impact of the relativistic effects in the cosmological
observables is therefore more subtle than some power increase $\sim({\cal
H}/k)^{n>0}$ at low $k$ (see [30] for the infrared behaviors in Fourier space,
which are still non-trivial). A sufficient set of conditions for this low
$k$-power cancellation are diffeomorphism invariance and adiabatic
fluctuations at large enough scales, both of which hold in the standard
cosmology, or the $\Lambda$CDM model. Our conclusions are supported by several
explicit perturbative computations at linear order [21, 31, 32, 33, 30, 34,
35], to which we recently added the first detailed non-linear example [36].
Our present demonstration, however, does not restrict to a particular
observable and applies to all orders in perturbation theory.

Infrared (in)sensitivity of the cosmological observables in linear theory.— In
the Newtonian theory the cosmological observables only depend on the
gravitational potential $\phi$ through its second spatial derivatives
$\partial_{i}\partial_{j}\phi\propto\partial_{i}\partial_{j}\mathcal{R}$ ,
i.e. the matter density contrast $\delta_{m}$ (trace part) and the tidal
tensor (traceless part). In the relativistic theory, however, some
contributions to the cosmological observables are directly proportional to the
metric fluctuations $h_{\mu\nu}\propto\mathcal{R}$ (e.g., the Sachs-Wolfe
effect [37]), including gravitational waves $\gamma_{ij}$, so they contribute
with lower powers of $k$ in the kernels of the observable $N$-point statistics
[14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 38, 29].

At the linear order in perturbations, the only non-trivial relation for the
observable $N$-point statistics is Eq. (2). In the Newtonian theory,
$\langle\partial^{2}\phi\,\partial^{2}\phi\rangle\sim
k^{4}P_{\mathcal{R}}(k)$, so $\mathcal{T}_{l}\sim{\cal O}(k^{2})$, while the
relativistic effects contribute at $\mathcal{T}_{l}\sim{\cal O}(k^{0\leq
n<2})$, and $\langle\phi\phi\rangle\sim P_{\mathcal{R}}(k)$ at worst.
Incidentally, the local primordial non-Gaussianity in galaxy clustering also
induces such low-derivative terms through the bias model [39, 40, 41], i.e., a
term $\langle\delta_{m}\phi\rangle\sim k^{2}P_{\mathcal{R}}(k)$. Now since
$P_{\mathcal{R}}(k)\sim k^{n_{s}-4}$ at low $k$, the most infrared-sensitive
${\cal O}(k^{0})$ contributions in $\mathcal{T}_{l}$ lead to an infrared-
divergent integral in Eq. (2) for non-blue spectra $n_{s}\leq 1$, as
observations currently favor (see, e.g., [42, 43]). At the level of the
angular power spectrum $C_{l}(z,z^{\prime})$, this divergence is relatively
harmless, because it arises only in the monopole $C_{0}$ [17, 18, 19, 20, 21,
38] (and also in the quadrupole $C_{2}$ if one includes tensor modes
$\gamma_{ij}$ [44]), which is usually ignored in literature. However, since
all multipoles contribute to the correlation function in Eq. (1), the presence
of infrared divergences in any multipoles leads to inconsistencies in the
observable two-point correlation function. Moreover, it was shown [45, 46, 47]
that ignoring the monopole fluctuation $C_{0}$ in the CMB or supernova
observations because of such divergences biases the cosmological parameter
estimation and underestimates their uncertainties.

This issue was resolved in a series of recent works for the luminosity
distance [31, 32], galaxy clustering [21, 33, 30, 35] and cosmic microwave
background [34], which showed that with a fully gauge-invariant expression for
the cosmological observables the infrared-sensitive contributions actually
cancel each other out. More precisely, by assuming linear Gaussian
$\Lambda$CDM and comoving observer/sources with the free-falling matter fluid,
and by expanding the kernel $\mathcal{T}_{l}(k)$ in $k$ one finds that all
terms $\propto k^{0\leq n<2}$ sum up to zero. Thus, the cosmological
observable $O(z,\hat{n})$ in the relativistic theory is just as infrared-
insensitive as its Newtonian counterpart $\mathcal{T}_{l}\sim{\cal O}(k^{2})$.
Importantly, for this cancellation to occur the fluctuations must be adiabatic
in the infrared limit [34, 35], which is the case in $\Lambda$CDM. Indeed,
setting $\lim_{k\to 0}k^{n}\mathcal{T}^{2}_{l}(k)=0$ for $n\in\\{0,1,2,3\\}$
one finds the equations satisfied by Weinberg’s adiabatic mode [48], i.e., the
$k\to 0$ limit of adiabatic solutions.

Another important condition for this cancellation to take place is that one
must consistently include all the relativistic contributions to the observable
$O(z,\hat{n})$, so that the coordinate-independence of the theoretical
expression is satisfied [21, 31, 32, 49, 33, 30, 34, 35]. The general form of
a cosmological observable at the linear order in perturbations is

$O(z,\hat{n})=\bar{O}(z)\left[1+X_{\rm o}(\hat{n})+X_{\rm
los}(z,\hat{n})+X_{\rm s}(z,\hat{n})\right]\,,$ (3)

where $\bar{O}$ is the background contribution, $X_{\rm o}$ represents
quantities at the observer point, $X_{\rm los}$ field fluctuations integrated
over the background line-of-sight path and $X_{\rm s}$ fields evaluated at the
source point 222Note that this separation is somewhat arbitrary, as one can
move some terms around by integrating by parts the line-of-sight terms.
Moreover, although each individual contribution can be gauge-dependent, the
sum $O(z,\hat{n})$ is invariant under any coordinate transformation at each
order, i.e., both the usual gauge transformations of cosmological perturbation
theory and more general elements of the diffeomorphism group [54].. In the
standard practice the “observer terms” $X_{\rm o}$ are simply discarded (see
below), thus yielding infrared-sensitive results. However, since $X_{\rm o}$
has no position dependence $(z,\hat{n})$, it can only contribute up to the
first few multipoles $C_{s}$, where $s=0,1,2$ is the highest spin of the
fields/quantities appearing in $X_{\rm o}$, which is why the spurious infrared
divergences only occur for the first few multipoles

Infrared (in)sensitivity in non-linear theory.— At the non-linear level the
field combinations appearing in the integrals of $O(z,\hat{n})$ generalize as
follows: the dominant ones at small scales are of the form $\sim
h^{n}\partial^{2}h$, while the more infrared-sensitive (lower-derivative) ones
are $\sim h^{n}$, $h^{n}\partial h$ and $h^{n}\partial h\partial h$, i.e. all
possible combinations with up to two spatial derivatives (from light
propagation and field equations). The first example of infrared insensitivity
at non-linear order was shown for the observed galaxy number density
bispectrum [36]: in the squeezed equal-redshift configuration $\langle
O(z,\hat{n})\,O(z,\hat{n})\,O(z,-\hat{n})\rangle$ the leading-order
relativistic effect contributions in the small momentum variable $k_{l}$
cancel out. The computation is long and non-trivial, but the busy reader can
still have a look at Eq. (4.32) of [36], where each contribution to the signal
is labelled by its physical interpretation and discussed. One can then see
explicitly the cancellation between light-propagation effects and in
particular the necessity of including the observer terms for this to happen.

An important aspect of non-linear perturbation theory is that, unlike the
linear case, the relativistic effects at the observer position affect all
multipoles of all observable multi-spectra. Indeed, at non-linear order there
are also cross-type terms in $O(z,\hat{n})$, e.g., of the form $X_{\rm
o}\times X_{\rm s}(z,\hat{n})\subset O^{(2)}(z,\hat{n})$, thus contributing
terms like $\langle[X_{\rm o}X_{\rm s}(z,\hat{n})]\,[X_{\rm o}X_{\rm
s}(z^{\prime},\hat{n}^{\prime})]\rangle$ in the 2-point function and
$\langle[X_{\rm o}X_{\rm s}(z,\hat{n})]\,X_{\rm
s}(z^{\prime},\hat{n}^{\prime})\,X_{\rm
s}(z^{\prime\prime},\hat{n}^{\prime\prime})\rangle$ in the 3-point function.
There is a priori no reason why the corresponding contributions would be
smaller than the ones without any observer dependence, especially at large
scales. The observer term issue is therefore not at all trivial non-linearly,
yet ignored in most existing computations of the bispectrum. Although the non-
linear realm currently lacks a more general rigorous exploration, i.e.
including all terms at a given order, one can still expect cancellations of
the lowest powers of the wave-number dependencies that are expected from
individual relativistic effects. Indeed, as we will show below, the underlying
reason is the diffeomorphism symmetry, which is independent of the
perturbative order.

General proof of infrared insensitivity.— Let us now demonstrate the infrared
insensitivity of relativistic effects in the cosmological observables. We
start by noting that cosmological observations involve a finite patch of
space-time – the one containing the observer-source system – whose typical
length-scale we denote by $L_{z}$. Let us then choose a reference scale $L\gg
L_{z}$ and split the metric perturbation $h_{\mu\nu}$ into long and short
wavelength contributions $h_{\mu\nu}\equiv h_{\mu\nu}^{\rm L}+h_{\mu\nu}^{\rm
S}$, in some given coordinates. One possibility for this split is the Gaussian
coarse-graining

$h_{\mu\nu}^{\rm L}(t,x):=\frac{1}{(2\pi)^{3/2}L^{3}}\int{\rm
d}^{3}y\exp\left[-\frac{|x-y|^{2}}{2L^{2}}\right]h_{\mu\nu}(t,y)\,,$ (4)

or in Fourier space,

$h_{\mu\nu}^{\rm
L}(t,k)\equiv\exp\left[-\frac{1}{2}\,L^{2}k^{2}\right]h_{\mu\nu}(t,k)\,.$ (5)

As a result, the short-mode metric fluctuation is

$h_{\mu\nu}^{\rm S}(t,k)\equiv h_{\mu\nu}(t,k)-h_{\mu\nu}^{\rm
L}(t,k)\sim{\cal O}(k^{2}h_{\mu\nu})\,,$ (6)

so the infrared-sensitive ${\cal O}(k^{0\leq n<2})$ contributions to the
cosmological observable $O(z,\hat{n})$ would disappear if $h_{\mu\nu}^{\rm
L}=0$.

By construction $h^{\rm L}_{\mu\nu}$ varies very little within a distance
$L_{z}$ from the observer world-line $x^{i}_{\rm o}(t)$, since we chose $L\gg
L_{z}$. We can therefore approximate the long-mode metric fluctuation by its
first-order expansion around the world-line

$h^{\rm L}_{\mu\nu}(t,\bm{x})\approx f_{\mu\nu}(t)+\left[x^{i}-x^{i}_{\rm
o}(t)\right]f_{i,\mu\nu}(t)\,.$ (7)

Note that the corrections ${\cal O}(\Delta x^{2}/L^{2})$ can be made
arbitrarily small by increasing $L$ and thus negligible for the purpose of
computing $O(z,\hat{n})$, as it only involves fields within a distance $L_{z}$
from the world-line. Put differently, using a large, yet finite $L\gg L_{z}$,
the approximated long-mode metric in Eq. (7) yields an observable that is
indistinguishable from the exact one. Incidentally, note that the second-order
derivative terms $\sim h^{n}\partial^{2}h$ are precisely those that vanish in
Eq. (7), while the infrared-sensitive ones ($h^{n}$, $h^{n}\partial h$ and
$h^{n}\partial h\partial h$) do not.

The central observation now is that, under some conditions discussed below,
the precise form in Eq. (7) can be removed by a residual large coordinate
transformation [48, 51, 52, 53], i.e., one that preserves the chosen gauge and
therefore the expression of the observable in terms of the perturbation
fields. Since cosmological observable relations are fully coordinate-
independent [54, 55], they are invariant under this manipulation of a
coordinate transformation, so the result of their computation must be the same
when $h_{\mu\nu}^{\rm L}=0$. Given that this is the part of $h_{\mu\nu}$ that
carries the infrared-sensitive terms ${\cal O}(k^{0\leq n<2})$ in the
integrals of the observable $O(z,\hat{n})$, we conclude that the corresponding
infrared-sensitive contributions for the observable $N$-point statistic must
be absent 333Note that one can always eliminate $h^{\rm L}_{\mu\nu}$ in Eq.
(7), i.e. for all values of $f_{\mu\nu}$ and $f_{i,\mu\nu}$, by a coordinate
transformation such as the one to conformal Fermi coordinates [60], but the
resulting expression of the observable would be different in terms of the
perturbations in the new coordinates, hence the need for a residual gauge
transformation in our argument..

In the typical gauges chosen in the literature one has $h_{\mu\nu}^{\rm L}\neq
0$, so the absence of low $k$-powers can only occur through cancellations
among infrared-sensitive terms. Since the argument above relies on the
diffeomorphism symmetry, one must include all the relativistic contributions
in the cosmological observables at a given order to maintain its coordinate-
independence, and this is not possible if some contributions – for example at
the observer position – are ignored. Intuitively, including those terms at the
observer position in $O(z,\hat{n})$ can then be understood as consistently
implementing the finite observer-source separation $L_{z}$ that sets a
reference scale with respect to which infrared modes are filtered.

Required conditions.— The long-mode part given by Eq. (7) can be eliminated
through a residual coordinate transformation

$h_{\mu\nu}\to\tilde{h}_{\mu\nu}\,,\hskip 28.45274pt\tilde{h}_{\mu\nu}^{\rm
L}=0\,,$ (8)

only for specific profiles $h_{\mu\nu}^{\rm L}$, i.e. specific functions
$f_{\mu\nu}(t)$ and $f_{i,\mu\nu}(t)$. Solutions that become
“pure-(residual)-gauge” in the infrared limit are known as “adiabatic modes”
[48] (see [57] for an up-to-date systematic exploration). Equivalently,
adiabatic modes can be generated by acting on purely short-mode solutions with
specific diffeomorphisms that preserve the gauge. At the linear level, where
all modes decouple, this means that they can be generated by acting with some
coordinate transformation on the background solution and in the matter sector
this leads to the relations

$\frac{\delta\rho^{a}}{\dot{\bar{\rho}}^{a}}=\frac{\delta
p^{a}}{\dot{\bar{p}}^{a}}=\frac{\delta\rho^{\rm tot}}{\dot{\bar{\rho}}^{\rm
tot}}=\frac{\delta p^{\rm tot}}{\dot{\bar{p}}^{\rm tot}}\,,$ (9)

which is thus consistent with the usual definition of “adiabatic modes” in
linear cosmological perturbation theory. Given Eq. (6), an adiabatic
$h_{\mu\nu}^{\rm L}$ means that any non-adiabatic contribution in $h_{\mu\nu}$
must decay at least as fast as $\sim k^{2}$ in the infrared.

An important subtlety is that acting with a large diffeomorphism can also
generate unphysical solutions. The typical example in linearized general
relativity (GR) [48] is the constant dilatation $x^{i}\to(1+\lambda)x^{i}$ in
the Newtonian gauge which transforms the potentials as $\phi\to\phi+\lambda$
and $\psi\to\psi$. The resulting solution no longer satisfies $\phi=\psi$ at
$k=0$ (assuming no anisotropic stress), which is allowed since that equation
is actually $\partial^{2}\phi=\partial^{2}\psi$. Adiabatic modes therefore
include the constraint of being smooth $k\to 0$ limits of finite-$k$ solutions
and, as a result, their profile depends on the equations of motion of the
theory. For instance, if we consider a modified gravity theory, then there
would still exist adiabatic modes, but the corresponding functions
$f_{\mu\nu}(t)$ and $f_{i,\mu\nu}(t)$ of Eq. (7) could be different from those
in GR, since the equations of motion would be different. However, most
modified gravity theories contain GR as a subset of solutions (e.g., Horndeski
theory solutions with constant scalar), so in that case one would at least
have the adiabatic modes of GR.

Let us now mention the most relevant adiabatic modes for GR in the context of
cosmology. There is the already mentioned scalar Weinberg adiabatic mode [48,
53], which is constant in space, but also its tensor and conformal group
generalization that includes a linear $\sim x^{i}$ term [51, 52, 53]. This is
exactly the spatial dependence we need to reproduce the metric in Eq. (7) for
cosmological applications. More precisely, it has been shown that one can
generate the Newtonian-gauge scalar components [48, 53] and both scalar and
tensor components in the comoving-gauge [51, 52]. There are two subtleties in
the presence of tensor modes: the linear $\sim x^{i}$ scalar profile can be
generated only if there are no short-mode tensor modes $\tilde{\gamma}_{ij}=0$
[51, 52] and one can only generate $\gamma_{ij}^{\rm L}$ to linear order in
that field [52].

The condition that only the adiabatic scalar degree of freedom is active above
some scale is satisfied in $\Lambda$CDM, since this property is satisfied in
the initial conditions and preserved through evolution non-linearly at scales
above the largest sound horizon among the involved propagating fields. Thus,
choosing $L$ above the largest sound horizon and $L_{z}$, we conclude that in
$\Lambda$CDM all cosmological observables are infrared-insensitive to all
orders in perturbation theory.

The condition of adiabaticity at the linear level was already pointed out in
[34] through its standard form in Eq. (9) and in [35] through the appearance
of the equations describing Weinberg’s adiabatic mode. Intuitively, deviations
from adiabaticity correspond to more than one physical reference frame, such
as the case of two distinct bulk flows influencing the observer-source system.
Even if one eliminates their average effect (adiabatic mode) through the
equivalence principle, their coordinate-independent difference (entropy mode)
will remain (see, e.g., [58, 59]). The complications that arise in the
presence of tensor modes, as described above, could also be understood as a
deviation from adiabaticity, i.e. an extra physical reference that obstructs
or remains after the elimination of a constant force.

Finally, so far we have implicitly assumed that all relativistic effects are
expressible in terms of $h_{\mu\nu}$, which is the case only if one assumes
free-falling observer and sources. In the explicit computations [21, 31, 32,
33, 30, 34, 35] both were taken to follow the matter fluid, itself being in
free-fall, which is accurate at the large scales of interest. However, this
free-fall assumption for the observer and sources is not necessary, thanks to
the requirement of adiabaticity at large enough scales. Independently of
whether these actors are subject to non-gravitational forces, or follow
different non-free-falling fluids, the important point is that their dynamics
at cosmological scales are ultimately determined by some of the available
fields in the theory of interest. These fields will then enter in the
expression for the observable through the observer/source perturbations.
However, since only the single (adiabatic) degree of freedom survives at the
largest scales, the infrared part of the extra fields will be expressible in
terms of $h_{\mu\nu}$ as well.

The conceptual issue with observer terms— We argued that consistent statistics
require the inclusion of observer terms in $O(z,\hat{n})$. However, for their
contribution to take the required form, i.e. proportional to primordial
spectra (2), observer terms must be expressed as stochastic fields evaluated
at the observer position. This is achieved by relating the observer quantities
to the underlying fields that govern the dynamics (e.g. $\bm{v}_{\rm
o}=\bm{v}_{m}(t_{\rm o},\bm{x}_{\rm o})$ for a comoving observer). Thus, the
observer adapts to each universe realization in the ensemble average, as the
sources do.

That approach, however, raised objections in private communications and
workshops, which revealed more concretely why observer terms were ignored for
so long in the literature. Critics would point out that the observer terms
$X_{\rm o}$ should not be considered as stochastic, as some of them can in
principle be measured with local experiments, in contrast to the rest of the
light-cone information, which is available only in a model-dependent way.

Importantly, there can be no “right” and “wrong” approach, since these are
just two different prescriptions for doing statistics with our single
universe. Which part of the information we model as deterministic or
stochastic is a matter of choice, similarly to the question of prior knowledge
in Bayesian statistics. Thus, the fact that some terms in $X_{\rm o}$ could be
measured directly in principle does not force one to fix that value. The
question therefore is: which prescription is the most convenient, most
effective, or even simply computable?

The local measurements required in the fixed-observer prescription are
dominated by the highly non-linear local dynamics, which one does not resolve
within perturbation theory, so their inclusion in this context would be
ambiguous. Moreover, although some quantities at the observer are measurable
(e.g. our relative velocity to the CMB), most of them are not (e.g. the scalar
potentials), as they cannot be expressed in terms coordinate/frame-independent
quantities (e.g. local curvature invariants). In contrast, the “statistical-
observer” prescription advocated here requires no extra input and is
straightforwardly computable. The condition of a comoving observer (and
sources) with matter is accurate at the scales one resolves within
perturbation theory and it is also an unambiguous, physical condition.
Moreover, the consistency of this approach was demonstrated in [38], where we
showed that, if we proceed without privileging the observer position at any
step, then the variance of unbiased estimators is exactly cosmic variance,
i.e. the usual information bound in cosmology.

Discussion.— We saw that infrared insensitivity is guaranteed if non-adiabatic
fluctuations decay at least as fast as $\sim k^{2}$ at low $k$. This condition
is intriguing, because one could a priori think of instances where it does not
hold, yet breaking it too much would lead to infrared-divergent observables –
clearly an unphysical result for a computation involving a finite space-time
patch. For instance, this happens at the linear level if non-adiabatic
contributions do not decay ($\mathcal{T}_{l}\sim k^{0}$) and $n_{s}\leq 1$.
One is therefore led to wonder whether there exists any mechanism that can
produce non-adiabatic fluctuations that survive at $k=0$. Lacking any definite
answer on this matter, we can at least say that this does seem at odds with
the underlying framework of cosmological perturbation theory, as one would
expect homogeneity and isotropy at $k\to 0$, thus leaving only room for
solutions that are coordinate artefacts in that limit, i.e. adiabatic
solutions.

###### Acknowledgements.

We are grateful to the participants of the workshop “General Relativistic
effects in observing the Large Scale Structure of the Universe” (June 2022,
Porto) for useful discussions. JY acknowledges useful discussions with Robert
Brandenberger. EM is supported by a Forschungskredit Grant of the University
of Zurich (grant FK-21-130) and JY is supported by the Swiss National Science
Foundation and a Consolidator Grant of the European Research Council.

## References

  * Peebles [1980] P. J. E. Peebles, _The large-scale structure of the universe_ (Princeton University Press, Princeton, 1980).
  * Peacock [1999] J. A. Peacock, _Cosmological Physics_ (1999).
  * Dodelson [2003] S. Dodelson, _Modern cosmology_ (2003).
  * Weinberg [2008] S. Weinberg, _Cosmology_ (Oxford University Press, ISBN 0-198-52682-2, 2008).
  * Note [1] The Fourier coordinates $k^{i}$ considered here should not be confused with the ones observers use to map the sky, as the two types are equivalent only in the Newtonian limit. Indeed, the “observer’s Fourier space” is the dual of Cartesian coordinates inferred from the data $(z,\hat{n})$ through some fiducial cosmology and thus parametrizing the past light-cone, not a constant-time hypersurface. Consequently, in this observational Fourier space large scales correspond to early times, meaning that there is either no well-defined infrared limit (finite particle horizon), or that limit is trivial since structure is absent in the distant past.
  * Percival et al. [2001] W. J. Percival, C. M. Baugh, J. Bland-Hawthorn, T. Bridges, R. Cannon, S. Cole, M. Colless, C. Collins, W. Couch, G. Dalton, et al., Mon. Not. R. Astron. Soc. 327, 1297 (2001), eprint astro-ph/0105252.
  * Eisenstein et al. [2005] D. J. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, et al., Astrophys. J. 633, 560 (2005), eprint arXiv:astro-ph/0501171.
  * Stubbs et al. [2004] C. W. Stubbs, D. Sweeney, J. A. Tyson, and LSST Collaboration, in _American Astronomical Society Meeting Abstracts_ (2004), vol. 36 of _Bulletin of the American Astronomical Society_ , p. 108.02.
  * Dewdney et al. [2009] P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. W. Lazio, IEEE Proceedings 97, 1482 (2009).
  * Laureijs et al. [2011] R. Laureijs, J. Amiaux, S. Arduini, J. . Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, et al., ArXiv e-prints (2011), eprint 1110.3193.
  * Green et al. [2012] J. Green, P. Schechter, C. Baltay, R. Bean, D. Bennett, R. Brown, C. Conselice, M. Donahue, et al. (2012), eprint 1208.4012.
  * Levi et al. [2013] M. Levi, C. Bebek, T. Beers, R. Blum, R. Cahn, D. Eisenstein, B. Flaugher, K. Honscheid, R. Kron, O. Lahav, et al., ArXiv e-prints (2013), eprint 1308.0847.
  * Doré et al. [2014] O. Doré, J. Bock, M. Ashby, P. Capak, A. Cooray, R. de Putter, T. Eifler, N. Flagey, Y. Gong, S. Habib, et al., arXiv e-prints arXiv:1412.4872 (2014), eprint 1412.4872.
  * Sasaki [1987] M. Sasaki, Mon. Not. R. Astron. Soc. 228, 653 (1987).
  * Futamase and Sasaki [1989] T. Futamase and M. Sasaki, Phys. Rev. D 40, 2502 (1989).
  * Bonvin et al. [2006] C. Bonvin, R. Durrer, and M. A. Gasparini, Phys. Rev. D 73, 023523 (2006), eprint arXiv:0511183.
  * Yoo et al. [2009] J. Yoo, A. L. Fitzpatrick, and M. Zaldarriaga, Phys. Rev. D 80, 083514 (2009), eprint arXiv:0907.0707.
  * Yoo [2010] J. Yoo, Phys. Rev. D 82, 083508 (2010), eprint arXiv:1009.3021.
  * Challinor and Lewis [2011] A. Challinor and A. Lewis, Phys. Rev. D 84, 043516 (2011), eprint arXiv:1105.5292.
  * Bonvin and Durrer [2011] C. Bonvin and R. Durrer, Phys. Rev. D 84, 063505 (2011), eprint arXiv:1105.5280.
  * Jeong et al. [2012] D. Jeong, F. Schmidt, and C. M. Hirata, Phys. Rev. D 85, 023504 (2012), eprint arXiv:1107.5427.
  * Ben-Dayan et al. [2014] I. Ben-Dayan, R. Durrer, G. Marozzi, and D. J. Schwarz, Phys. Rev. Lett. 112, 221301 (2014), eprint 1401.7973.
  * Umeh et al. [2014] O. Umeh, C. Clarkson, and R. Maartens, Classical and Quantum Gravity 31, 202001 (2014), eprint 1207.2109.
  * Yoo and Zaldarriaga [2014] J. Yoo and M. Zaldarriaga, Phys. Rev. D 90, 023513 (2014), eprint 1406.4140.
  * Bertacca et al. [2014] D. Bertacca, R. Maartens, and C. Clarkson, J. Cosmol. Astropart. Phys. 9, 037 (2014), eprint 1405.4403.
  * Di Dio et al. [2014] E. Di Dio, R. Durrer, G. Marozzi, and F. Montanari, J. Cosmol. Astropart. Phys. 12, 017 (2014), eprint 1407.0376.
  * Bonvin et al. [2015] C. Bonvin, C. Clarkson, R. Durrer, R. Maartens, and O. Umeh, J. Cosmol. Astropart. Phys. 6, 50 (2015), eprint 1503.07831.
  * Koyama et al. [2018] K. Koyama, O. Umeh, R. Maartens, and D. Bertacca, J. Cosmol. Astropart. Phys. 2018, 050 (2018), eprint 1805.09189.
  * Magi and Yoo [2022] M. Magi and J. Yoo, J. Cosmol. Astropart. Phys. 2022, 071 (2022), eprint 2204.01751.
  * Grimm et al. [2020] N. Grimm, F. Scaccabarozzi, J. Yoo, S. G. Biern, and J.-O. Gong, J. Cosmol. Astropart. Phys. 2020, 064 (2020), eprint 2005.06484.
  * Biern and Yoo [2017a] S. G. Biern and J. Yoo, J. Cosmol. Astropart. Phys. 4, 045 (2017a), eprint 1606.01910.
  * Biern and Yoo [2017b] S. G. Biern and J. Yoo, J. Cosmol. Astropart. Phys. 026 (2017b), eprint 1704.07380.
  * Scaccabarozzi et al. [2018] F. Scaccabarozzi, J. Yoo, and S. G. Biern, J. Cosmol. Astropart. Phys. 10, 024 (2018), eprint 1807.09796.
  * Baumgartner and Yoo [2021] S. Baumgartner and J. Yoo, Phys. Rev. D 103, 063516 (2021), eprint 2012.03968.
  * Castorina and Di Dio [2022] E. Castorina and E. Di Dio, J. Cosmol. Astropart. Phys. 2022, 061 (2022), eprint 2106.08857.
  * Yoo et al. [2022] J. Yoo, N. Grimm, and E. Mitsou, J. Cosmol. Astropart. Phys. 2022, 050 (2022), eprint 2204.03002.
  * Sachs and Wolfe [1967] R. K. Sachs and A. M. Wolfe, Astrophys. J. 147, 73+ (1967).
  * Mitsou et al. [2020] E. Mitsou, J. Yoo, R. Durrer, F. Scaccabarozzi, and V. Tansella, Physical Review Research 2, 033004 (2020), eprint 1905.01293.
  * Dalal et al. [2008] N. Dalal, O. Doré, D. Huterer, and A. Shirokov, Phys. Rev. D 77, 123514 (2008), eprint 0710.4560.
  * Matarrese and Verde [2008] S. Matarrese and L. Verde, Astrophys. J. Lett. 677, L77 (2008), eprint 0801.4826.
  * Slosar et al. [2008] A. Slosar, C. Hirata, U. Seljak, S. Ho, and N. Padmanabhan, J. Cosmol. Astropart. Phys. 8, 31 (2008), eprint 0805.3580.
  * Planck Collaboration et al. [2016] Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, et al., Astron. Astrophys. 594, A13 (2016), eprint 1502.01589.
  * Planck Collaboration et al. [2020] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, et al., Astron. Astrophys. 641, A6 (2020), eprint 1807.06209.
  * Adamek et al. [2016] J. Adamek, R. Durrer, and V. Tansella, J. Cosmol. Astropart. Phys. 1, 024 (2016), eprint 1510.01566.
  * Yoo et al. [2019a] J. Yoo, E. Mitsou, N. Grimm, R. Durrer, and A. Refregier, J. Cosmol. Astropart. Phys. 2019, 015 (2019a), eprint 1905.08262.
  * Yoo et al. [2019b] J. Yoo, E. Mitsou, Y. Dirian, and R. Durrer, Phys. Rev. D 100, 063510 (2019b), eprint 1905.09288.
  * Yoo [2020] J. Yoo, Phys. Rev. D 101, 043507 (2020), eprint 1911.07869.
  * Weinberg [2003] S. Weinberg, Phys. Rev. D 67, 123504 (2003), eprint astro-ph/0302326.
  * Yoo and Durrer [2017] J. Yoo and R. Durrer, J. Cosmol. Astropart. Phys. 9, 016 (2017), eprint 1705.05839.
  * Note [2] Note that this separation is somewhat arbitrary, as one can move some terms around by integrating by parts the line-of-sight terms. Moreover, although each individual contribution can be gauge-dependent, the sum $O(z,\hat{n})$ is invariant under any coordinate transformation at each order, i.e., both the usual gauge transformations of cosmological perturbation theory and more general elements of the diffeomorphism group [54].
  * Creminelli et al. [2012] P. Creminelli, J. Noreña, and M. Simonović, J. Cosmol. Astropart. Phys. 2012, 052 (2012), eprint 1203.4595.
  * Hinterbichler et al. [2012] K. Hinterbichler, L. Hui, and J. Khoury, J. Cosmol. Astropart. Phys. 2012, 017 (2012), eprint 1203.6351.
  * Creminelli et al. [2013] P. Creminelli, J. Noreña, M. Simonović, and F. Vernizzi, J. Cosmol. Astropart. Phys. 2013, 025 (2013), eprint 1309.3557.
  * Mitsou and Yoo [2020] E. Mitsou and J. Yoo, Springer Briefs in Physics arXiv:1908.10757 (2020), eprint 1908.10757.
  * Mitsou and Yoo [2022] E. Mitsou and J. Yoo, Classical and Quantum Gravity 39, 217003 (2022), eprint 2205.01544.
  * Note [3] Note that one can always eliminate $h^{\rm L}_{\mu\nu}$ in Eq. (7), i.e. for all values of $f_{\mu\nu}$ and $f_{i,\mu\nu}$, by a coordinate transformation such as the one to conformal Fermi coordinates [60], but the resulting expression of the observable would be different in terms of the perturbations in the new coordinates, hence the need for a residual gauge transformation in our argument.
  * Pajer and Jazayeri [2018] E. Pajer and S. Jazayeri, J. Cosmol. Astropart. Phys. 2018, 013 (2018), eprint 1710.02177.
  * Marozzi et al. [2013] G. Marozzi, G. P. Vacca, and R. H. Brandenberger, J. Cosmol. Astropart. Phys. 2, 027 (2013), eprint 1212.6029.
  * Brandenberger et al. [2018] R. Brandenberger, L. L. Graef, G. Marozzi, and G. P. Vacca, Phys. Rev. D 98, 103523 (2018), eprint 1807.07494.
  * Dai et al. [2015] L. Dai, E. Pajer, and F. Schmidt, J. Cosmol. Astropart. Phys. 2015, 043 (2015), eprint 1502.02011.