\section{Conclusions and discussion}
We have proposed a new variational inference framework for estimating the intensity of a continuous sigmoidal Cox process. By seeing an augmented input space from a superposition of two \acro{ppp}{s},  we have derived a scalable and computationally efficient structured variational approximation. Our framework does not require discretization or accurate numerical computation of integrals on the input space, it is not limited to specific kernel functions and properly accounts for the strong dependencies existing across the latent variables.
Through extensive empirical evaluation we have shown that our methods compares favorably  against `exact'  but computationally costly \acro{mcmc} schemes, while being almost an order of magnitude faster. More importantly, our inference scheme outperforms all competing approaches in terms of uncertainty quantification. The benefit of the proposed scheme and resulting \acro{svi} are particularity pronounced on multivariate input settings where accounting for the highly coupled variables become crucial for interpolation and prediction.
Future work will focus on relaxing the factorization assumption between the \acro{gp} and the latent points. Introducing a fully structured variational inference would further improve the accuracy performance of the method but would require further approximations in the variational objective.

\section{Structured Variational Inference in the augmented space}
Given the joint distribution in Eq.\xspace \ref{eq:joint_lik}, our goal is to estimate the posterior distribution over all latent variables given the data. i.e.\xspace $p(\vec{f},\vec{u},M, \{\vec{y}_m\}_{m=1}^{M},\lambda^\star|\mathcal{D})$. This posterior is analytically intractable and we resort to variational inference \cite{jordan1999introduction}. Variational  inference entails defining an approximate posterior $q(\vec{f},\vec{u},M, \{\vec{y}_m\}_{m=1}^{M},\lambda^\star)$  and optimizing the so-called evidence lower bound (\acro{elbo}) with respect to this distribution.
In \acro{sgcp}, the \acro{gp} and the latent variables are highly coupled and breaking their dependencies would lead to poor approximations, especially in high dimensional settings. Fig.\xspace \ref{fig:plates} shows the structure of a general posterior distribution for \acro{sgcp} without any factorisation assumption. We consider an approximate posterior distribution that takes dependencies into account:
\begin{align}
Q(\vec{f}, \vec{u}, M, \{\vec{y}_m\}_{m=1}^{M},\lambda^\star)  = p(\vec{f}|\vec{u}) q(\{\vec{y}_m\}_{m=1}^{M}|M) q(M|\vec{u}, \lambda^\star) q(\vec{u})q(\lambda^\star)
\label{eq:variational_posterior}
\end{align}

With respect to the general posterior distribution, the only factorisation we impose in Eq.\xspace \eqref{eq:variational_posterior} is in the factor $q(\{\vec{y}_m\}_{m=1}^{M}|M)$ where we drop the dependency on $\vec{f}$, see dashed line in Fig.\xspace \ref{fig:plates}. We set:
\begin{align*}
	& q(\vec{u}) = \mathcal{N}(\vec{m}, \mat{S}) \quad  \quad q(\lambda^\star) = \text{Gamma}(\alpha, \beta) \quad q(\{\vec{y}_m\}_{m=1}^{M}|M) = \prod_{m=1}^{M}\sum_{s=1}^{S}\pi_s \mathcal{N}_T(\mu_s, \sigma^2_s; \mathcal{X})
\end{align*}
 where $\mathcal{N}_T(\cdot;\mathcal{X})$ denotes a truncated Gaussian distribution on $\mathcal{X}$.
More importantly, we assume $q(M|\vec{u}, \lambda^\star) = \text{Poisson}(\eta)$ with $\eta = \lambda^\star\int_{\mathcal{X}}\sigma(-\vec{u}(\vec{x}))d\vec{x}$. This is indeed the \emph{true} conditional posterior distribution for the number of thinned points, see Proposition (3.7) in \cite{moller2003statistical}. Considering $q(M|\vec{u}, \lambda^\star)$ we thus fully account for the dependency structure existing among $M$, $\vec{u}$ and $\lambda^\star$.
Crucially, our algorithm does not require \emph{accurate} estimation of $\int_{\mathcal{X}}\sigma(-\vec{u}(\vec{x}))d\vec{x}$.  Differently from the competing techniques, where the algorithm convergence and the posterior $q(\vec{f})$ is \emph{directly} dependent on numerical integration, we only require evaluation of the integral during the optimisation but $q(\vec{f})$ and thus $\lambda(\vec{x})$ do not directly depend on its value.

\subsection{Evidence Lower Bound}
Following standard variational inference arguments, it is straightforward to show that the \acro{elbo} decomposes as:
\begin{align*}
\calL_{\text{elbo}}
& = N (\psi(\alpha) - \log(\beta)) - V\frac{\alpha}{\beta} - \log(N!) +  \underbrace{\expectation{Q}{M \log(\lambda^\star)}}_{T_1}   -  \underbrace{\expectation{Q}{\log(M!)}}_{T_2} +
\\& +\sum_{n=1}^{N} \expectation{q(\vec{f})}{\log(\sigma(\vec{f}(\vec{x}_n)))}
+ \underbrace{\expectation{Q} {\sum_{m=1}^{M} \log(\sigma(-\vec{f}(\vec{y}_m)))}}_{T_3} - \calL_{\text{kl}}^{\vec{u}}- \calL_{\text{kl}}^{\lambda^\star} -  \underbrace{\calL_{\text{ent}}^M}_{T_4}- \underbrace{\calL_{\text{ent}}^{\{\vec{y}_m\}_{m=1}^M}}_{T_5} \numberthis
\label{eq:elbo}
\end{align*}
where $V = \int_{\mathcal{X}}d\vec{x}$, $\psi(\cdot)$ is the digamma function and $q(\vec{f}) = \mathcal{N}(\mat{A}\vec{m},\mat{K}_{xx} - \mat{A}\mat{K}_{zx} +\mat{A} \mat{S} \mat{A}')$ . The terms denoted by $T_i, i=1,...,5$ cannot be computed analytically. Na\"{i}vely, black-box variational inference algorithms could be used to estimate these terms via Monte Carlo, thus sampling from the full variational posterior (Eq.\xspace \eqref{eq:variational_posterior}). This would require sampling $\vec{f}$, $\lambda^\star$, $M$ and $\{\vec{y}_m\}_{m=1}^M$ thus slowing down the algorithm while leading to slow convergence. On the contrary, we exploit the structure of the model and the approximate posterior to simplify these terms and increase the algorithm efficiency. Denote $\mu(\u) = \int_{\mathcal{X}}{\sigma(-\vec{u}(\vec{x}))d\vec{x}}$, we can write:

\begin{align}
	&T_1 = \expectation{q(\lambda^\star)}{\lambda^\star\log(\lambda^\star)} \expectation{q(\vec{u})}{\mu(\u)}
	\quad  T_2 \leq \left[\expectation{q(\lambda^\star)}{\lambda^\star\log(\lambda^\star)} + \frac{\alpha}{\beta}\right]\expectation{q(\vec{u})}{\mu(\u)\log(\mu(\u))} \label{eq:derivations1} \\
	&T_3  = \frac{\alpha}{\beta}\expectation{q(\vec{u})}{\mu(\u)}\expectation{q(\vec{f})q(\vec{y}_m)}{\log(\sigma(-\vec{f}(\vec{y}_m)))} \quad \quad T_5 =\frac{\alpha}{\beta} \expectation{q(\vec{y}_m)}{\log q(\vec{y}_m)} \expectation{q(\vec{u})}{\mu(\u)} \label{eq:derivations2}\\
	&T_4 \leq \frac{\alpha}{\beta}\expectation{q(\vec{u})}{\mu(\u)(\log(\mu(\u))-1)} + \expectation{q(\lambda^\star)q(\vec{u})}{\lambda^\star\log(\lambda^\star)\expectation{q(\vec{u})}{\mu(\u)-1}}
	\label{eq:derivations}
\end{align}
where we use Stirling's approximation for $\log M! \approx M\log(M) - M$ and the bound $\eta \log(\eta) \leq \expectation{}{M\log(M)} \leq \eta \log(\eta + 1)$. See the supplement (\S1) for the full derivations.  Eqs.\xspace \eqref{eq:derivations1}--\eqref{eq:derivations} give a lower bound for $\calL_{\text{elbo}}$ which avoids sampling from $ q(M|\vec{u}, \lambda^\star)$ and $q(\{\vec{y}_m\}_{m=1}^{M}|M)$ and does not require computing the \acro{gp} on the stochastic locations. The remaining expectations are with respect to reparameterizable distributions. We thus avoid the use of score function estimators which would lead to high-variance gradient estimates. Stochastic optimisation techniques can be used to evaluate $T_3$ and $\sum_{n=1}^{N} \expectation{q(\vec{f})}{\log(\sigma(\vec{f}(\vec{x}_n)))}$ making the corresponding computational cost independent of $M$ and $N$.  This would reduce the computational complexity of the algorithm to $\mathcal{O}(K^3)$. However, when the number of inputs used per mini-batch equals $N$, the time complexity becomes $\mathcal{O}(NK^2)$. In the following experiments, we show how the proposed structured approach together with these efficient \acro{elbo} computations lead to higher predictive performances and better uncertainty quantification.

\section{Introduction}

Point processes have been used effectively to model a variety of event data such as occurrences of diseases \cite{diggle2013spatial, lloyd-icml-2015}, location of earthquakes \cite{marsan2008extending} or crime events \cite{grubesic2008spatio, aglietti2018} . The most commonly adopted class of models for such discrete data are non-homogenous Poisson processes and in particular Cox processes \cite{cox1955some}. In these, the observed events are assumed to be generated from a Poisson point process (\acro{ppp}) whose intensity is stochastic, enabling non-parametric inference and uncertainty quantification.

Gaussian processes \citep[\acro{gp}{s};][]{Rasmussen:2005:GPM:1162254} form a flexible prior over functions and, therefore, have been used to model the intensity of a Cox process via a non-linear positive link function. Typical mappings are the exponential \cite{moller1998log, diggle2013spatial}, the square \cite{lloyd-icml-2015,pmlr-v37-lian15} and the sigmoidal \cite{adams2009tractable, gunter2014efficient, donner2018efficient} transformations.  In general, inferring the intensity function over a continuous input space $\mathcal{X}$ is highly problematic and different algorithms have been proposed to deal with this issue depending on the transformation used. For example, under the exponential transformation, a regular computational grid is commonly introduced \cite{diggle2013spatial}. While this significantly simplifies inference, it leads to poor approximations, especially in high dimensional settings. Increasing the resolution of the grid to improve the approximation  yields computationally prohibitive algorithms that do not scale, highlighting the well-known trade-off between statistical performance and computational cost.

Other algorithms have been proposed to deal with a continuous $\mathcal{X}$ but they are computationally expensive \cite{adams2009tractable, gunter2014efficient}, are limited to simple covariance functions \cite{lloyd-icml-2015}, require accurate numerical integration over the domain \cite{donner2018efficient} or do not account for the model dependencies in the posterior distribution \cite{donner2018efficient}.  In this paper we propose an inference framework that addresses all of these modeling and inference limitations by having a tractable representation of the likelihood via augmentation with a superposition of \acro{ppp}{s}. This enables a scalable structured variational inference algorithm (\acro{svi}) in the continuous space directly, where the approximate posterior distribution incorporates dependencies between the variables of interest. Our specific contributions are as follows.

\textbf{Scalable inference in continuous input spaces}: The augmentation of the input space via a process superposition view allows us to develop a scalable variational inference algorithm that does not require discretization or \emph{accurate} numerical integration. With this view, we obtain a joint distribution that is readily normalized, providing a natural regularization over the latent variables in our model.

\textbf{Efficient structured posterior estimation}: We estimate a joint posterior  that captures the complex variable dependencies in the model while being significantly faster than sampling approaches.

\textbf{State-of-the-art performance}: Our experimental evaluation shows the benefits of our approach when compared to state-of-the-art  inference schemes, link functions, augmentation schemes and representations of the input space $\mathcal{X}$.

\section{Model formulation}
We consider learning problems where we are given a dataset of $N$ events  $\mathcal{D} = \{\vec{x}_n\}_{n=1}^{N}$, where $\vec{x}_n$  is a $d$-dimensional vector in the compact space $\mathcal{X} \subset \mathbb{R}^D$. We aim at modeling these data via a \acro{ppp}, inferring the latent intensity function $\lambda(\vec{x}):\mathcal{X} \to \mathbb{R}^+$ and making probabilistic predictions.

\subsection{Sigmoidal Gaussian Cox process}
Consider a realization $\xi = (N, \{\vec{x}_1,...,\vec{x}_n\})$ of a \acro{ppp} on $\mathcal{X}$ where the points $\{\vec{x}_1, ..., \vec{x}_n\}$ are treated as \textit{indistinguishable} apart from their locations \cite{daley2007introduction}.
Conditioned on $\lambda(\vec{x})$, the Cox process likelihood function evaluated at $\xi$ can be written as:
\begin{align}
	\mathcal{L}(\xi|\lambda(\vec{x})) = \exp\left(-\int_{\tau}\lambda(\vec{x})d\vec{x}\right) \prod_{n=1}^{N}\lambda(\vec{x}_n) ,
\label{eq:intractable_likelihood}
\end{align}
where the intensity is given by $\lambda(\vec{x}) = \lambda^\star\sigma(f(\vec{x}))$ with $\lambda^\star >0$ being an upperbound on $\lambda(\vec{x})$,  $\sigma(\cdot)$ denoting the the logistic sigmoid function and $f$ is drawn from a zero-mean \acro{gp} prior with covariance function $\covfunc{\vec{x}}{\x^{\prime}}{\vecS{\theta}}$ and hyperparameters $\vecS{\theta}$, i.e.\xspace $f | \vecS{\theta} \sim \mathcal{GP}(\vec{0}, \covfunc{\vec{x}}{\x^{\prime}}{\vecS{\theta}})$.  We will refer to this joint model as the sigmoidal Gaussian Cox process (\acro{sgcp} ).  Notice that, when considering the tuple $(\vec{x}_1, ..., \vec{x}_n)$ instead of the set $\{\vec{x}_1, ..., \vec{x}_n\}$, and thus the event $\xi_0= (N, (\vec{x}_1,...,\vec{x}_n))$, the likelihood function is given by  $\mathcal{L}(\xi_0|\lambda(\vec{x}))  = \frac{\mathcal{L}(\xi|\lambda(\vec{x}))}{N!}$. There are indeed $N!$ permutations of the events $\{\vec{x}_1,...,\vec{x}_n\}$ giving the same point process realization. When the set $\{\vec{x}_1, ..., \vec{x}_n\}$ is known, considering $\mathcal{L}(\xi|\lambda(\vec{x}))$ or $\mathcal{L}(\xi_0|\lambda(\vec{x}))$ does not affect the inference procedure. The same holds for \acro{mcmc} algorithms inferring the event locations. In this case, the factorial term disappears in the computation of the acceptance ratio. However, as we shall see later, when the event locations are latent variables in a model and inference proceeds via a variational approximation the difference between the two likelihoods is essential. Indeed, while $\mathcal{L}(\xi_0|\lambda(\vec{x}))$ is normalized with respect to $N$, one must be cautious when integrating the likelihood in Eq.\xspace \ref{eq:intractable_likelihood} over sets and bring back the missing $N!$ factor so as  to obtain a proper discrete probability mass function for $N$.

As it turns out, inference in \acro{sgcp} is \textit{doubly intractable}, as it requires solving the  integral in Eq.\xspace \eqref{eq:intractable_likelihood} and then computing the intractable  posterior distribution $p(\vec{f}_N, \lambda^\star| \{\vec{x}_n\}_{n=1}^{N})$, where $\vec{f}_N$ denotes the latent function at the $N$ event locations. One way to achieve tractability is through augmentation of the input space \cite{adams2009tractable,donner2018efficient}, a procedure that introduces precisely those latent (event) variables that require explicit normalization during variational inference.  We will describe below a process superposition view of this augmented scheme that allows us to define a proper distribution over the joint space of observed and latent variables and carry out posterior estimation via variational inference.

\subsection{Augmentation via superposition}
A very useful property of independent \acro{ppp}{s} is that their superposition, which is defined as the combination of events from two processes in a single one, is a \acro{ppp}. Consider two \acro{ppp} with intensities $\lambda(\vec{x})$ and $\nu(\vec{x})$ and realisations $(N, \{\vec{x}_1, ...,\vec{x}_n \})$ and $(M, \{\vec{y}_1, ..., \vec{y}_M\})$
respectively. The combined event $\xi_R = (R = M+N, \{\vec{v}_1, ..., \vec{v}_{R}\})$ is a realization of a \acro{ppp} with intensity given by $\lambda(\vec{x}) + \nu(\vec{x}))$ where knowledge of which points originated from which process is assumed lost. The likelihood for $\mathcal{L}(\xi_R|\lambda(\vec{x}),\nu(\vec{x}))$ can be thus written as:
\begin{equation}
	\sum_{N=0}^{R} \binom{N + M}{N} \sum_{P_N \in \mathbb{P}_N}\left(\frac{\exp(-\int_{\mathcal{X}} \lambda(\vec{x})d\vec{x})}{N!}\prod_{r \in P_N}\lambda(r)\times\frac{\exp(-\int_{\mathcal{X}} \nu(\vec{x})d\vec{x})}{M!}\prod_{r \in P_N^c}\nu(r)\right),
\end{equation}
where $\mathbb{P}_N$ denotes the collection of all possible partitions of size $N$, $P_N$ represents an element of $\mathbb{P}_N$ and $P_N^c$ is its complement.

Consider now $R=N+M$ to be the total number of events resulting from thinning \cite{lewis1979simulation} where $N$ is the number of observed events while $M$ is the number of latent events with stochastic locations $\vec{y}_1, ..., \vec{y}_M$. We assume that the probability of observing an event is given by $\sigma(f(\vec{x}))$ while the probability for the event to be latent is $\sigma(-f(\vec{x}))$. In addition, let $\lambda^\star\int_{\mathcal{X}}d\vec{x}$ be the expected total number of events. We can see the realization $(M+N, (\vec{x}_1, ..., \vec{x}_N, \vec{y}_1, ..., \vec{y}_M))$ as the result of the superposition of two \acro{ppp}{s} with intensities $\lambda(x) = \lambda^\star\sigma(f(\vec{x}))$ and $\nu(\vec{x}) = \lambda^\star\sigma(-f(\vec{x}))$. Differently from the standard superposition, we do know which events are observed and which are latent. In writing the likelihood for $(M+N, \{\vec{x}_1, ..., \vec{x}_N, \vec{y}_1, ..., \vec{y}_M\})$ we thus do not need to consider all the possible partitions of $N$. We can write $\mathcal{L}_{N+M}   \stackrel{\text{\tiny def}}{=}  \mathcal{L}(N+M, (x_1, ..., x_N, y_1, ..., y_M))$:
\begin{align}
\mathcal{L}_{N+M}
	&= \frac{\exp(-\int_{\mathcal{X}} \lambda(\vec{x})d\vec{x})}{N!}\prod_{r \in P_N}\lambda(r)\times\frac{\exp(-\int_{\mathcal{X}} \nu(\vec{x})d\vec{x})}{N!}\prod_{r \in P_N^c}\nu(r) \\
	&= \frac{1}{N!M!}\exp(-\lambda^\star\int_{\mathcal{X}}dx)(\lambda^\star)^{M+N} \prod_{n=1}^{N} \sigma(\vec{f}(\vec{x}_n))\prod_{m=1}^{M} \sigma(-\vec{f}(\vec{x}_m)).
	\label{eq:tuples_extended_lik}
\end{align}
There is a crucial difference between Eq.\xspace \ref{eq:tuples_extended_lik} and the usual likelihood considered in \acro{sgcp}. Eq.\xspace \ref{eq:tuples_extended_lik} represents a distribution over tuples and thus, as mentioned above, is properly normalized. In addition, it makes a distinction between the observed and latent events and it is thus different from Eq.\xspace \ref{eq:intractable_likelihood} written for the the tuple $(M+N, \{\vec{x}_1, ..., \vec{x}_N, \vec{y}_1, ..., \vec{y}_M\})$.
We can write the full joint distribution as $\mathcal{L}^{+}_{N+M} \stackrel{\text{\tiny def}}{=} \mathcal{L}(\{\vec{x}_n\}_{n=1}^{N}, \{\vec{y}_m\}_{m=1}^{M}, M, \vec{f}, \lambda^\star| \tau, \vecS{\theta})$:
\begin{equation}
\label{eq:joint_lik}
 \mathcal{L}^{+}_{N+M} = \frac{(\lambda^\star)^{N+M}\exp(-\lambda^\star\int_{\mathcal{X}}dx) }{N!M!}
\prod_{n=1}^{N} \sigma(\vec{f}(\vec{x}_n))  \prod_{m=1}^{M} \sigma(-\vec{f}(\vec{y}_m)) \times p(\vec{f}) \times p(\lambda^\star)  ,
\end{equation}
where $p(\vec{f}) \stackrel{\text{\tiny def}}{=}  p(\vec{f}_{N+M}) $ denotes the joint prior at both  $\{\vec{x}_n\}_{n=1}^N$ and  $\{\vec{y}_m\}_{m=1}^M$ and $p(\lambda^\star)$ denotes the prior over the upper bound of the intensity function. We consider $p(\lambda^\star) = \text{Gamma}(a, b)$ and set $a$ and $b$ so that $\lambda^\star$ as has mean and standard deviation equal to 2$\times$ and 1$\times$ the intensity we would expect  from an homogenous Poisson process on $\mathcal{X}$.

\subsection{Scalability via inducing variables}
As in standard \acro{gp} modulated models, the introduction of a \acro{gp} prior poses significant computational challenges during posterior estimation as inference would be dominated by algebraic operations that are cubic on the number of observations.  In order to make inference scalable, we follow the inducing-variable approach proposed
by \cite{2009variational} and further developed by \cite{bonilla2016generic}. To this end, we consider an augmented prior $p(\vec{f}, \vec{u})$ with $K$ underlying inducing variables denoted by $\vec{u}$.  The corresponding inducing inputs are given by the $K \times D$ matrix $\mat{Z}$. Major computational gains are realized when $K \ll N+M$.
The augmented prior distributions for the inducing variables and the latent functions are $p(\vec{u} | \vecS{\theta}) = \mathcal{N}(\vec{0}, \mat{K}_{zz}) $ and $p( \vec{f}| \vec{u}, \vecS{\theta}) = \mathcal{N}(\mat{K}_{xz}  (\mat{K}_{zz})^{-1} \vec{u},  \mat{K}_{xx} - \mat{A} \mat{K}_{zx}) $ where $\mat{A} = \mat{K}_{xz} (\mat{K}_{zz})^{-1}$.
The  matrices $\mat{K}_{xx}$, $\mat{K}_{xz}$, $\mat{K}_{zx}$ and $(\mat{K}_{zz})^{-1}$ are the covariances induced by evaluating the corresponding covariance functions at all pairwise rows of the event locations $\{\vec{x}_n,\vec{y}_m\}_{n=1,m=1}^{N, M}$ and the inducing inputs $\mat{Z}$.
\input{plate}

\subsection{Related work}\label{sec:related_work}
\acro{gp}-modulated point processes are the gold standard for modeling event data. Performing inference in these models, e.g.\xspace under the exponential transformation, has typically required discretization where the domain is gridded and the intensity function is assumed to be constant over each grid cell \cite{moller1998log, cunningham2008fast, diggle2013spatial,brix2001spatiotemporal}. Alternatively, \citet{lasko2014efficient} also considers an exponential link function and performs inference over a renewal process resorting to numerical integration  within a computationally expensive sampling scheme.
These methods  suffer from poor scaling with the dimensionality of $\mathcal{X}$ and sensitivity to the choice of discretization or numerical integration technique. Several approaches have been proposed to deal with inference in the continuous domain directly by using alternative transformations along with additional modeling assumptions and computational tricks or by constraining the \acro{gp} \cite{lopez2019gaussian}.

One of those alternative transformation is the squared mapping as developed in the Permanental process \citep{pmlr-v70-walder17a,lloyd-icml-2015,pmlr-v37-lian15, lloyd2016latent, john2018large}. Although the square transformation enables analytical computation of the required integrals over $\mathcal{X}$, this only holds for certain standard types of kernels such as the squared exponential. In addition, Permanental processes suffer from important identifiability issues such as reflection invariance and lead to model with ``nodal lines'' \cite{john2018large}.

Another transformation is the scaled logistic sigmoid function proposed by \cite{adams2009tractable}, which exploits augmentation of the input space via thinning \cite{lewis1979simulation} to achieve tractability. This model is known as the sigmoidal Gaussian Cox process (\acro{sgcp}). Their proposed inference algorithm is based on Markov chain Monte Carlo (\acro{mcmc}), which enables drawing `exact' samples from the posterior intensity. However, as acknowledged by the authors,  it has significant computational demands making it inapplicable to large datasets. As an extension to this work,  \cite{gunter2014efficient} introduce the concept of ``adaptive thinning'' and propose an expensive \acro{mcmc} scheme which scales as $\mathcal{O}(N^3)$.
More recently, \cite{donner2018efficient} introduced a neat double augmentation scheme for \acro{sgcp} which enables closed form updates using a mean-field approximation (\acro{vi}-\acro{mf}).  However, it requires accurate numerical integration over $\mathcal{X}$, which makes the performance of the algorithm highly dependent on the number of integration points.

In this work, we overcome the limitations of the mentioned \acro{vi}-\acro{mf} and \acro{mcmc} schemes by proposing an \acro{svi} framework, henceforth \acro{stvb}, which takes into account the complex posterior dependencies while being scalable and thus applicable to high-dimensional real-world settings.
To the best of our knowledge we are the first to propose a fast structured variational inference framework for \acro{gp} modulated point process models.
See Tab.\xspace \ref{tab:summary_related} for a summary of the most relevant related works.

\input{related_table}

\section{ELBO derivations}
Here we derive the expressions given in Eqs.\xspace (7)-(10). As given in Eq.\xspace (7) the evidence lower bound ($\calL_{\text{elbo}}$) decomposes as:

\begin{align*}
& \calL_{\text{elbo}} = \expectation{Q}{\log\left[\frac{p(\{\vec{x}_n\}_{n=1}^{N}, \{\vec{x}_m\}_{m=1}^{M}, M, \vec{f}, \vec{u}, \lambda^\star| \tau, \vecS{\theta})}{p(\vec{f}|\vec{u}) q(\{\vec{x}_m\}_{m=1}^{M}|M) q(M|\vec{u}, \lambda^\star) q(\vec{u})q(\lambda^\star)}\right]}
\\& = \expectation{Q}{\log{p(\{\vec{x}_n\}_{n=1}^{N}, \{\vec{x}_m\}_{m=1}^{M}, M, \vec{u},\lambda^\star| \tau, \vecS{\theta})}}
- \expectation{Q}{\log q(\{\vec{x}_m\}_{m=1}^{M}|M) q(M|\vec{u}, \lambda^\star) q(\vec{u})q(\lambda^\star)}
\\& =\mathbb{E}_Q((N + M)\log(\lambda^\star) - \lambda^\star \mu(\tau) - \log(M!) - \log(N!) + \sum_{n=1}^{N}\log(\sigma(\vec{f}(\vec{x}_n))) \\& +  \sum_{m=1}^{M}\log(\sigma(-\vec{f}(\vec{x}_m)) +  \log(p(\vec{u}))   +  \log(p(\lambda^\star))
\\ & - \log(q(\vec{u})) - \log(q(M|\vec{u}, \lambda^\star)) - \log(q(\lambda^\star)) - \log(q(\{\vec{x}_m\}_{m=1}^{M}|M))
\\& = N (\psi(\alpha)-\log(\beta)) - V\frac{\alpha}{\beta} - \log(\log N!) + \underbrace{\expectation{Q}{M \log(\lambda^\star)}}_{T_1}  -  \underbrace{\expectation{Q}{\log M!}}_{T_2}  +
\\& +\sum_{n=1}^{N} \expectation{q(\vec{u})}{\log(\sigma(\vec{f}(\vec{x}_n)))}
+ \underbrace{\expectation{Q} {\sum_{m=1}^{M} \log(\sigma(-\vec{f}(\vec{x}_m)))}}_{T_3}+
\\& - KL(q(\vec{u})||p(\vec{u})) - KL(q(\lambda^\star)||p(\lambda^\star))
\\& -  \underbrace{\expectation{Q}{\log q(M|\vec{u}, \lambda^\star)}}_{T_4} -  \underbrace{\expectation{Q}{\log q(\{x_m\}_{m=1}^{M}|M)}}_{T_5}
\end{align*}

Let's now focus on the terms $T_i, i=1, ..., 5$. In the following derivations we will exploit the Stirling's approximation for $\log(M!)\approx M\log M - M$ and the bounds $\eta \log \eta \leq \expectation{Q}{M \log M}\leq \eta \log(\eta +1)$ to obtain a lower bound for the $\calL_{\text{elbo}}$. The term $T_1$ (Eq.\xspace (8)) is given by:
\begin{align*}
T_1 &= \expectation{q(\vec{u})q(\lambda^\star)}{\expectation{q(M|\vec{u}, \lambda^\star)}{M\log(\lambda^\star)}}
\\& = \expectation{q(\vec{u})q(\lambda^\star)}{\log(\lambda^\star)  \expectation{q(M|\vec{u}, \lambda^\star)}{M}}
\\& = \expectation{q(\vec{u})q(\lambda^\star)}{\log(\lambda^\star)  \lambda^\star \int_{\tau}\sigma(-\vec{u}(x))dx}
\\& = \expectation{q(\lambda^\star)}{\lambda^\star\log(\lambda^\star)} \expectation{q(\vec{u})}{\mu(\u)}
\end{align*}

The term $T_2$ (Eq.\xspace (8)) is given by:

\begin{align*}
T_2 & \approx \expectation{Q}{M\log(M)} - \expectation{Q}{M}
\\& = \expectation{Q}{M\log(M)} - \expectation{q(\lambda^\star)q(\vec{u})}{\expectation{q(M|\lambda^\star, \vec{u})}{M}}
\\& = \expectation{Q}{M\log(M)} - \expectation{q(\lambda^\star)q(\vec{u})}{\lambda^\star  \int_{\tau}\sigma(-\vec{u}(x))dx}
\\& = \expectation{Q}{M\log(M)} - \expectation{q(\lambda^\star)}{\lambda^\star } \expectation{q(\vec{u})}{\int_{\tau}\sigma(-\vec{u}(x))dx}
\\& = \expectation{Q}{M\log(M)} - \frac{\alpha}{\beta} \expectation{q(\vec{u})}{\mu(\u)}
\\& \geq \expectation{q(\lambda^\star)q(\vec{u})}{\eta \log (\eta +1)} - \frac{\alpha}{\beta} \expectation{q(\vec{u})}{\mu(\u)}
\\& = \expectation{q(\lambda^\star)}{\lambda^\star\log(\lambda^\star)}\expectation{q(\vec{u})}{\mu(\u)} + \frac{\alpha}{\beta}\expectation{q(\vec{u})}{\mu(\u)\log(\mu(\u))} - \frac{\alpha}{\beta} \expectation{q(\vec{u})}{\mu(\u)}
\end{align*}

The term $T_3$ (Eq.\xspace (9)) is given by:
\begin{align*}
T_3 &= \expectation{q(\vec{f})q(\vec{u})q(\vec{y}_m)q(\lambda^\star)}{\expectation{q(M|\vec{f}, \vec{y}_m)}{\sum_{m=1}^{M}\log(\sigma(-\vec{f}(\vec{y}_m)}}
\\& =  \expectation{q(\vec{f})q(\vec{u})q(\vec{y}_m)q(\lambda^\star)}{\log(\sigma(-\vec{f}(\vec{y}_m) \expectation{q(M|\vec{f}, \vec{y}_m)}{\sum_{m=1}^{M}1}}
\\& =  \expectation{q(\vec{f})q(\vec{u})q(\vec{y}_m)}{\log(\sigma(-\vec{f}(\vec{y}_m) \lambda^\star \mu(\u)}
\\& =\frac{\alpha}{\beta}\expectation{q(\vec{u})}{\mu(\u)}\expectation{q(\vec{f})q(\vec{y}_m)}{\log(\sigma(-\vec{f}(\vec{y}_m)))}
\end{align*}

The term $T_4$ (Eq.\xspace (10)) is given by:

\begin{align*}
T_4 &= \expectation{Q}{-\lambda^\star\mu(\u)} + \expectation{Q}{M\log(\lambda^\star\mu(\u)}  -\expectation{Q}{\log M!}
\\& = - \frac{\alpha}{\beta}\expectation{q(\vec{u})}{\int_{\tau}\sigma(-\vec{u}(x))dx}
+ \expectation{q(\lambda^\star)q(\vec{u})}{\lambda^\star\log(\lambda^\star)\mu(\u) + \lambda^\star\mu(\u)\log(\mu(\u))}\expectation{q(\vec{u})}{\int_{\tau}\sigma(-\vec{u}(x))dx}
\\& - \expectation{Q}{\log M!}
\\& = - \frac{\alpha}{\beta}\expectation{q(\vec{u})}{\int_{\tau}\sigma(-\vec{u}(x))dx} +\expectation{q(\lambda^\star)}{\lambda^\star\log(\lambda^\star)}\expectation{q(\vec{u})}{\mu(\u)}
\\& + \frac{\alpha}{\beta}\expectation{q(\vec{u})}{\mu(\u)\log(\mu(\u))}
\\& + \frac{\alpha}{\beta} \expectation{q(\vec{u})}{\mu(\u)} - \expectation{Q}{M\log(M)}
\end{align*}

Substituting the bound for $\expectation{Q}{M\log(M)}$ and simplifying we get to Eq.\xspace (10).

Finally, the term $T_5$ (Eq.\xspace (9)) is given by:
\begin{align*}
T_5 & = \expectation{Q}{\sum_{m=1}^M \log q(\vec{y}_m)}
\\& = \expectation{q(\vec{u}, \lambda^\star, M)}{\sum_{m=1}^M \expectation{q(\vec{y}_m)}{\log q(\vec{y}_m)}}
\\& = \expectation{q(\vec{u}, \lambda^\star, M)}{M} \expectation{q(\vec{y}_m)}{\log q(\vec{y}_m)}
\\& = \expectation{q(\vec{u}) q(\lambda^\star)}{\lambda^\star \mu(\u)}\expectation{q(\vec{y}_m)}{\log q(\vec{y}_m)}
\\& = \frac{\alpha}{\beta}\expectation{q(\vec{y}_m)}{\log q(\vec{y}_m)} \expectation{q(\vec{u})}{\mu(\u)}
\end{align*}

\section{Performance metrics}

We test the algorithms evaluating the $\mathit{l}_2$ norm to the true intensity function (in the synthetic settings), the test log likelihood ($\ell_{test}$) on the test set and the negative log predicted likelihood (\acro{nlpl}) on the training set. These metrics are computed as follow:

\begin{align}
	\mathit{l}_2 = \int_{\mathcal{X}}(\lambda(\vec{x}) - \bar{\lambda}(\vec{x}))^2d\vec{x}
\end{align}
where $\lambda(\vec{x})$ is the true intensity function, $ \bar{\lambda}(\vec{x})$ is the posterior mean intensity and the integral is evaluated numerically.

\begin{align}
	\ell_{test} = \expectation{q(\lambda^\star)q(\vec{f})}{\log \left[\exp\left(-\int_{\mathcal{X}}\lambda(\vec{x})d\vec{x}\right)\prod_{\vec{x} \in \mathcal{D_{\text{test}}}}\lambda(\vec{x})\right]}
\end{align}

where again the integral is computed via numerical integration.

The \acro{nlpl} is computed as
\begin{align}
	\acro{nlpl} =  - \frac{1}{S}\sum_{s=1}^{S}\log p(N_{\text{train}}|\int_{\mathcal{X}}\lambda^s(\vec{x})d(\vec{x}))
\end{align}

where $S$ denotes the number of samples from the variational distributions $q(\vec{f})$ and $q(\lambda^\star)$.

Finally, the \acro{ec} is computed by evaluating the coverage of the \acro{ci}{s} of the posterior ($p(N | \mathcal{D})$) and predictive ($p(N^{*} | \mathcal{D})$). To construct the empirical count distribution we sample from the variational distributions $q(\vec{f})$ and $q(\mat{W})$, obtain samples of $\lambda(\vec{x})$ and simulate $N$ or $N^*$ from $\text{Poisson}(\lambda^\star  \int_{\mathcal{X}}\sigma(\vec{f}(\vec{x}))d\vec{x})$.

\section{Additional experimental results}
For all comparisons we consider a \acro{gp} with squared-exponential covariance function with equally set hyperparameters. Denote by $\vecS{\theta}_i = (l, \sigma^2_f)$ the values of the hyperameters for the kernel function $K(\vec{x}, \vec{x}')$ on $\lambda_i(\vec{x})$ where $l$ indicates the lenghtscale. We set:
\begin{itemize}
	\item $\vecS{\theta}_1 = (10,1)$
	\item $\vecS{\theta}_2 = (0.25,1)$
	\item $\vecS{\theta}_3 = (15,1)$
\end{itemize}

For the real-world settings we have:
\begin{itemize}
	\item $\vecS{\theta}_{\text{neuronal data}}= (10,1)$
	\item $\vecS{\theta}_{\text{taxi data}}= (0.3,1)$
\end{itemize}

\subsection{Synthetic data experiments}
In Tab.\xspace \ref{tab:coverage_appendix_1}, \ref{tab:coverage_appendix_2} and \ref{tab:coverage_appendix_3} we report the values of \acro{ec} for different \acro{ci}{s} and on both the training and test set.

\begin{table}
	\caption{$\lambda_1(\vec{x})$ - \acro{ec} performance on training and test dataset. Higher values are better.}
	\begin{center}
\resizebox{\columnwidth}{!}{\begin{tabular}{lcccccccccccccccccc}
	\toprule
	&  \multicolumn{5}{c}{$\lambda_1(\vec{x})$ - In-sample \acro{ec}}    &  \multicolumn{5}{c}{$\lambda_1(\vec{x})$ - Out-of-sample \acro{ec}}   \\
	\cmidrule(lr){2-6}
	\cmidrule(lr){7-11}

	& 10\% \acro{ci}             & 20\%    \acro{ci}             & 30\% \acro{ci}               & 40\%  \acro{ci}              & 50\%        \acro{ci}        & 10\%      \acro{ci}              & 20\%        \acro{ci}            & 30\%       \acro{ci}             & 40\%     \acro{ci}              & 50\%     \acro{ci}               \\
	\cmidrule(lr){2-6}
	\cmidrule(lr){7-11}
	\acro{stvb}   & \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}             &  \textbf{1.00}             & 0.96                 & 0.88                 & \textbf{0.81}                 & \textbf{0.72}                 & \textbf{0.60}                  \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.24)               & (0.24)               & (0.23)               & (0.29)               & (0.29)               \\
	\acro{mfvb} &  \textbf{1.00}              &  \textbf{1.00}               &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              & 0.95                 & 0.80                  & 0.76                 & 0.61                 & 0.52                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.00)                & (0.00)                & (0.00)                & (0.00)                & (0.00)                \\
	\acro{vbpp} &  \textbf{1.00}              &  \textbf{1.00}               &  \textbf{1.00}              &  \textbf{1.00}              & 0.10              & \textbf{1.00}                  & \textbf{0.97}                 & 0.75                 & 0.41                 & 0.04                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.30)            & (0.00)                & (0.05)               & (0.21)               & (0.25)               & (0.09)               \\
	\acro{sgcp}  &  \textbf{1.00}              &  \textbf{1.00}               &  \textbf{1.00}              &  \textbf{1.00}              & 0.60              & 0.75                 & 0.60                  & 0.39                 & 0.27                 & 0.08                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.49)           & (0.29)               & (0.33)               & (0.28)               & (0.22)               & (0.12)               \\
	\acro{lgcp}   & 0.70              & 0.00               & 0.00              & 0.00              & 0.00              & 0.48                 & 0.22                 & 0.08                 & 0.03                 & 0.01                 \\
	& (0.46)           & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.00)                & (0.00)                & (0.00)                & (0.00)                & (0.00)                \\
	\bottomrule
\end{tabular}}
		\label{tab:coverage_appendix_1}
\end{center}
\end{table}

\begin{table}
	\caption{$\lambda_2(\vec{x})$ - \acro{ec} performance on training and test dataset. Higher values are better.}
	\begin{center}
		\resizebox{\columnwidth}{!}{\begin{tabular}{lcccccccccccccccccccccc}
	\toprule
	&  \multicolumn{5}{c}{$\lambda_2(\vec{x})$ - In-sample \acro{ec}}    &  \multicolumn{5}{c}{$\lambda_2(\vec{x})$ - Out-of-sample \acro{ec}}   \\
\cmidrule(lr){2-6}
\cmidrule(lr){7-11}

	& 10\% \acro{ci}             & 20\%    \acro{ci}             & 30\% \acro{ci}               & 40\%  \acro{ci}              & 50\%        \acro{ci}        & 10\%      \acro{ci}              & 20\%        \acro{ci}            & 30\%       \acro{ci}             & 40\%     \acro{ci}              & 50\%     \acro{ci}               \\
	\hline
	\acro{stvb}    & \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}             &  \textbf{1.00}             & \textbf{1.00}                 & \textbf{0.97}                 & \textbf{0.91}                 & \textbf{0.88}                 & \textbf{0.86}                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.00)                & (0.09)               & (0.24)               & (0.23)               & (0.22)               \\
	\acro{mfvb} & \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}             &  \textbf{1.00}              & 0.92                 & 0.92                 & 0.89                 & 0.84                 & 0.82                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.00)                & (0.00)                & (0.00)                & (0.00)                & (0.00)                \\
	\acro{vbpp} & \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}                    & 0.10              & 0.92                 & 0.86                 & 0.76                 & 0.45                 & 0.05                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.30)            & (0.24)               & (0.23)               & (0.26)               & (0.26)               & (0.05)               \\
	\acro{sgcp}  &\textbf{1.00}            & 0.90               & 0.70              & 0.40              & 0.30              & 0.90                  & 0.90                  & 0.64                 & 0.14                 & 0.00                  \\
	& (0.00)            & (0.30)             & (0.46)           & (0.49)           & (0.46)           & (0.00)                & (0.00)                & (0.09)               & (0.05)               & (0.00)                \\
	\acro{lgcp}   & 0.10              & 0.00               & 0.00              & 0.00              & 0.00              & 0.80                  & 0.22                 & 0.04                 & 0.00                  & 0.00                  \\
	& (0.30)            & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.24)               & (0.16)               & (0.08)               & (0.00)                & (0.00)                \\
	\hline
\end{tabular}}
\label{tab:coverage_appendix_2}
\end{center}
\end{table}

\begin{table}
	\caption{$\lambda_3(\vec{x})$ - \acro{ec} performance on training and test dataset. Higher values are better.}
	\begin{center}
		\resizebox{\columnwidth}{!}{\begin{tabular}{lcccccccccccccccccccccc}
	\toprule
	&  \multicolumn{5}{c}{$\lambda_3(\vec{x})$ - In-sample \acro{ec}}    &  \multicolumn{5}{c}{$\lambda_3(\vec{x})$ - Out-of-sample \acro{ec}}   \\
\cmidrule(lr){2-6}
\cmidrule(lr){7-11}

& 10\% \acro{ci}             & 20\%    \acro{ci}             & 30\% \acro{ci}               & 40\%  \acro{ci}              & 50\%        \acro{ci}        & 10\%      \acro{ci}              & 20\%        \acro{ci}            & 30\%       \acro{ci}             & 40\%     \acro{ci}              & 50\%     \acro{ci}               \\
	\hline
	\acro{stvb}    &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}             &  \textbf{1.00}                & \textbf{1.00}                  & \textbf{1.00}                 & 0.99                 & 0.97                 & 0.92                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.00)                & (0.00)                & (0.00)                & (0.00)                & (0.12)               \\
	\acro{mfvb} &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}             &  \textbf{1.00}             & \textbf{1.00}                  & \textbf{1.00}                  & 0.97                 & 0.91                 & 0.78                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.00)                & (0.00)                & (0.00)                & (0.00)                & (0.00)                \\
	\acro{vbpp} &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}                    & 0.10              & 0.97                 & 0.94                 & 0.83                 & 0.43                 & 0.03                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.30)            & (0.09)               & (0.15)               & (0.19)               & (0.14)               & (0.05)               \\
	\acro{sgcp}  & 0.80              & 0.70               & 0.50              & 0.40              & 0.00              & 0.82                 & 0.54                 & 0.49                 & 0.34                 & 0.02                 \\
	& (0.40)            & (0.46)            & (0.50)            & (0.49)           & (0.00)            & (0.12)               & (0.05)               & (0.03)               & (0.07)               & (0.04)               \\
	\acro{lgcp}   & \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}              &  \textbf{1.00}             &  \textbf{1.00}           & \textbf{1.00}                  & \textbf{1.00}                  & \textbf{1.00}                  & \textbf{1.00}                  & \textbf{0.95}                 \\
	& (0.00)            & (0.00)             & (0.00)            & (0.00)            & (0.00)            & (0.00)                & (0.00)                & (0.00)                & (0.00)                & (0.00)                \\
	\hline
\end{tabular}}
\label{tab:coverage_appendix_3}
\end{center}
\end{table}

\subsection{Real data experiments}

In Tab.\xspace \ref{tab:coverage_real_appendix} we report the values of \acro{ec} for different \acro{ci}{s} and on both the training and test set.

\begin{table}
	\caption{Real data. Values are given as In-sample - Out-of-sample \acro{ec}. Mean and standard errors (in parenthesis) are computed across different seeds.}
	\begin{center}
		\begin{tabular}{@{\extracolsep{4pt}}@{\kern\tabcolsep}lcccccccccccccccc}

			\toprule
			& \multicolumn{5}{c}{Neuronal data} \\
			\cmidrule(lr){2-6}
			&  10\%  \acro{ec}&  20\% \acro{ec}& 30\%  \acro{ec}&  40\% \acro{ec}& 50\%\acro{ec} \\
			\hline
			\multirow{2}{*}{\acro{stvb}}    &  \textbf{1.00}-\textbf{1.00} &  \textbf{1.00}-\textbf{1.00} &  \textbf{1.00}-\textbf{1.00} & \textbf{0.99}-\textbf{0.56}      &  \textbf{0.01}- 0.00  \\
			&  0.00 - 0.00 & 0.00 - 0.00 &  0.00-0.00& (0.10)-(0.50) &  (0.10)-0.00 \\
			\multirow{2}{*}{\acro{mfvb}}   &   \textbf{1.00}-\textbf{1.00} &  \textbf{1.00}-0.62 &  \textbf{1.00}-0.03 & 0.78-0.00      &  0.00 - 0.00  \\
			& 0.00 - 0.00 &  0.00-(0.49)&  0.00-(0.17) & (0.41)-0.00 &  0.00 - 0.00 \\
			\multirow{2}{*}{\acro{vbpp}}   &    \textbf{1.00}-0.53 &  \textbf{1.00}-0.00 &  \textbf{1.00}-0.00 & 0.83-0.00    & \textbf{0.01}-0.00    \\
			& 0.00-(0.50) &  0.00 - 0.00&  0.00 - 0.00& (0.38)-0.00  &  (0.10)-0.00 \\
			\midrule
			& \multicolumn{5}{c}{Taxi data}\\
			\cmidrule(lr){2-6}
			&  10\% \acro{ec} &  20\%\acro{ec} & 30\%  \acro{ec}&  40\% \acro{ec}& 50\% \acro{ec}\\
			\hline
			\multirow{2}{*}{\acro{stvb}}    &  \textbf{1.00}-1.00      &  \textbf{1.00}-\textbf{1.00} & 0.81-\textbf{0.37}     & 0.09-\textbf{0.01}     &          0.00-0.00 \\
			&   0.00-0.00   &  0.00-0.00 & (0.39)-(0.48) & (0.29)-(0.10) &       0.00-0.00 \\
			\multirow{2}{*}{\acro{mfvb}}   &   0.49-0.93   &  0.00-0.13  & 0.00-0.00       & 0.00-0.00        &          0.00-0.00  \\
			&  (0.50)-(0.26) &  0.00-(0.34) & 0.00-0.00      &0.00-0.00       &    0.00-0.00 \\
			\multirow{2}{*}{\acro{vbpp}}     & \textbf{1.00}-0.00      &  \textbf{1.00}-0.00 & \textbf{0.98}-0.00     & \textbf{0.48}-0.00     &          0.00-  0.00  \\
			&  0.00- (0.00)     &  0.00-(0.00) & (0.14)-(0.00)     & (0.50)-0.00   &          0.00-0.00 \\
			\bottomrule
		\end{tabular}
	\label{tab:coverage_real_appendix}
	\end{center}
\end{table}

\end{document}

\section{Experiments}
We test our algorithm on three 1\acro{d} synthetic data settings and on two 2\acro{d} real-world applications\footnote{Code and data for all the experiments will be provided.}.

\textbf{Baselines}
We compare against alternative inference schemes, different link functions and a different augmentation scheme. In terms of continuous models, we consider the sampling approach of \citet{adams2009tractable}(\acro{sgcp}), a Permanental Point process model  \cite{lloyd-icml-2015}(\acro{vbpp}) and a mean-field  approximation based on a  P\'{o}lya-Gamma augmentation \cite{donner2018efficient}(\acro{mfvb}). In addition, we compare against a discrete variational log Gaussian Cox process model \cite{nguyen2014automated}(\acro{lgcp}). Details are given in the supplement (\S3).

\textbf{Performance measures}
We test the algorithms evaluating the $\mathit{l}_2$ norm to the true intensity function (for the synthetic datasets), the test log likelihood ($\ell_{test}$) on the test set and the negative log predicted likelihood (\acro{nlpl}) on the training set. In order to assess the model capabilities in terms of uncertainty quantification, we compute the empirical coverage (\acro{ec}), i.e.\xspace the coverage of the empirical count distributions obtained by sampling from the posterior intensity. We do that for different credible intervals (\acro{ci})  on both the training (in-sample, $p(N|\mathcal{D})$) and test set (out-of-sample, $p(N^*|\mathcal{D})$). Details on the metrics are in the supplement (\S2).
For the synthetic data experiments, we run the algorithms with 10 training datasets each including a different \acro{ppp} realization sampled from the ground truth.  For each different training set, we then evaluate the performance on other 10 unseen realizations sampled again from the ground truth. We compute the mean and the standard deviation for the presented metrics averaging across the training and test sets. For the real data settings, we compute the \acro{nlpl} and in-sample \acro{ec} on the observed events. We then test the algorithm computing both $\ell_{test}$ and out-of-sample \acro{ec} on the held-out events. In order to compute the out-of-sample \acro{ec} we rescale the intensity function as $\lambda_{test}(\vec{x})=  \lambda_{train}(\vec{x}) - N_{train}/V + N_{test}/V$ with $V = \int_{\mathcal{X}}d\vec{x}$. We then sample from $\lambda_{test}(\vec{x})$ and generate the predicted count distributions for different seeds.
\input{figure_synthetic}
\textbf{Synthetic experiments}
We test our approach using the three toy example proposed by \cite{adams2009tractable}:
\begin{itemize}
	\item $\lambda_1(\vec{x}) = 2\exp(-1/15) + \exp(-[(x-15)/10]^2)$ \quad  $x \in [0,50]$,
	\item $\lambda_2(\vec{x}) = 5\text{sin}(x^2) + 6$ \quad  $x \in [0,5]$ and
	\item $\lambda_3(\vec{x})$ piecewise linear through $(0,20)$, $(25,3)$,$(50,1)$, $(75,2.5)$, $(100,3)$ \quad  $x \in [0,100]$.
\end{itemize}
For \acro{lgcp}, we discretize the input space considering a grid cell width of one for $\lambda_1(\vec{x})$ and $\lambda_3(\vec{x})$ and of 0.5 for  $\lambda_2(\vec{x})$.  For \acro{mfvb} we consider 1000 integration points. In terms of $q(\{\vec{y}_m\}_{m=1}^M|M)$,  we set $S=5$ but consistent results where found across different values of this parameter.
The results are given in Fig.\xspace \ref{fig:new_synt} and Tab.\xspace \ref{tab:synthetic_1}, where we see that all algorithms recover similar predicted mean intensities and give roughly comparable performances across all metrics.  Out of all 9 settings and metrics (top section of Tab.\xspace \ref{tab:synthetic_1})  our method (\acro{stvb}) outperforms competing methods on 3 cases and it is only second to \acro{sgcp} on 6 cases. However, the \textsc{cpu} time  of \acro{sgcp} is almost an order of magnitude larger than ours even in these simple low-dimensional problems, making that method inapplicable to large datasets. This confirms the benefits of having structured approximate posteriors within  a computationally efficient inference algorithm such as \acro{vi}. In terms of uncertainty quantification (bottom section of Tab.\xspace \ref{tab:synthetic_1}), our algorithm outperforms  all competing approaches for $\lambda_1(\vec{x})$ and  $\lambda_2(\vec{x})$.
\input{table_synthetic}