Unlocking Point Processes through Point Set Diffusion

## Abstract

Abstract Point processes model the distribution of random point sets in mathematical spaces, such as spatial and temporal domains, with applications in fields like seismology, neuroscience, and economics.
Existing statistical and machine learning models for point processes are predominantly constrained by their reliance on the characteristic intensity function, introducing an inherent trade-off between efficiency and flexibility.
In this paper, we introduce Point Set Diffusion , a diffusion-based latent variable model that can represent arbitrary point processes on general metric spaces without relying on the intensity function.
By directly learning to stochastically interpolate between noise and data point sets, our approach enables efficient, parallel sampling and flexible generation for complex conditional tasks defined on the metric space.
Experiments on synthetic and real-world datasets demonstrate that Point Set Diffusion achieves state-of-the-art performance in unconditional and conditional generation of spatial and spatiotemporal point processes while providing up to orders of magnitude faster sampling than autoregressive baselines.

## 1 Introduction

Figure: Figure 1: Illustration of Point Set Diffusion for earthquakes in Japan. The forward process stochastically interpolates between the original data point set $X_{0}$ and a noise point set $X_{T}$, progressively removing data points and adding noise points. To generate new samples from the data distribution, we approximate the reverse posterior $q(X_{t}|X_{0},X_{t+1})$ and add approximate data points and remove noise points.
Refer to caption: x1.png

Point processes describe the distribution of point sets in a mathematical space where the location and number of points are random.
On Euclidean spaces, point processes (e.g., spatial and/or temporal; SPP, STPP, TPP) have been widely used to model events and entities in space and time, such as earthquakes, neural activity, transactions, and social media posts.

Point processes can exhibit complex interactions between points, leading to correlations that are hard to capture effectively (Daley & Vere-Jones, ).
The distribution of points is typically characterized by a non-negative intensity function, representing the expected number of events in a bounded region of space (Daley et al., ).
A common approach to modeling point processes on general metric spaces is to parameterize an inhomogeneous intensity as a function of space.
However, this approach assumes independence between points, which limits its ability to capture complex interactions and hinders generalization across different point sets (Daley et al., ; Daley & Vere-Jones, ).

For ordered spaces like time (STPP, TPP), the predominant approach is to model the conditional intensity autoregressively, where each point is conditioned on the past, allowing for temporal causal dependencies, which can be conveniently captured by state-of-the-art machine learning models Shchur et al. ().
While this enables point interactions, these models rely on likelihood-based training and autoregressive sampling, which require integrating the intensity function over the entire space.
Ultimately, this limits possible models, as it either necessitates oversimplified parameterizations that restrict point dependencies and introduce smoothness (Ozaki, ; Ogata, ; Zhou & Yu, ), or approximations with amortized inference (Zhou et al., ), numerical (Chen et al., ), or Monte Carlo methods (Hong & Shelton, ).
Thus, capturing complex point dependencies and sampling from point processes, particularly on general metric spaces, remains an open and challenging problem.

Lüdke et al. () overcame the limitations of the conditional intensity function for temporal point processes by proposing Add-Thin, a diffusion model for TPPs based on the thinning and superposition property of TPPs directly modeling entire event sequences.
In this paper, we generalize this idea to point processes on general metric spaces and derive a diffusion-based latent variable model, Point Set Diffusion, that directly learns to model the stochastic interpolation between a data point set and samples from any noise point process (see Figure 1).
Furthermore, we show how to generate conditional samples with our unconditional Point Set Diffusion model to solve arbitrary conditioning tasks on general metric spaces.
Our experiments demonstrate that Point Set Diffusion achieves state-of-the-art results on conditional and unconditional tasks for SPPs and STPPs. Our contributions can be summarized as follows:

- •
We derive a diffusion-based latent variable model for point processes on general metric spaces, capturing the distribution of arbitrary point processes by learning stochastic interpolations between data and noise point sets.
- •
Our model supports efficient and parallel sampling of point sets and enables generation for arbitrary conditional tasks defined as binary masks on the metric space.
- •
We introduce a model-agnostic generative evaluation framework for point process models on Euclidean spaces.
- •
Our method achieves state-of-the-art results for conditional and unconditional generation of SPPs and STPPs while offering orders of magnitude faster sampling.

## 2 Background

### 2.1 Point Processes

A point process (Daley et al., ) is a stochastic process where realizations consist of finite sets of randomly located points in a mathematical space. More formally, let $(D,d)$ be a complete, separable metric space equipped with its Borel $\sigma$-algebra $\mathcal{B}$. A point process on $D$ is a mapping $X$ from a probability space $(\Omega,\mathcal{A},\mathcal{P})$ into $N^{lf}$, the set of all possible point configurations, such that for any bounded Borel set $A\subseteq D$, the number of points in $A$, denoted by $N(A)$, is a finite random variable.

Given a realization of the point process $X=\{{\bm{x}}_{i}\in D\}_{1\leq i\leq n}$, where $n$ is the number of points, the number of points in a region is expressed as the counting measure $N(A)=\sum_{i=1}^{n}{\mathbf{1}\{{\bm{x}}_{i}\in A\}}$. Here, we assume the point process is simple, i.e., almost surely $N(\{{\bm{x}}_{i}\})\leq 1$ for all ${\bm{x}}_{i}\in D$, meaning no two points coincide.
Point processes are commonly characterized by their intensity function, which is defined through the following random measure:

$$ $A\mapsto\mu(A)\coloneqq\mathbb{E}[N(A)]=\int_{A}\lambda({\bm{x}})\,\mathop{}\! \mathrm{d}{\bm{x}},$ (1) $$

where $\mu(A)$ represents the expected number of points in a region $A$.
Then, a point process is said to have intensity $\lambda$ if the measure $\mu$ above has a density $\lambda$ with respect to the Lebesgue measure $\mu(A)$.
Thus, the intensity function $\lambda({\bm{x}})$ gives the expected number of points per unit volume in a small region of the Borel set $A\subseteq D$.

As the points in a realization $X$ can exhibit complex correlations, the intensity function can be non-trivial to parameterize.
For an Euclidean space $\mathbb{R}$ we can specify the Papangelou intensity (Daley et al., ):

$$ $\lambda({\bm{x}})=\lim_{\delta\rightarrow 0}\frac{P\{N(B_{\delta}({\bm{x}}))=1 |C(N(\mathbb{R}\setminus B_{\delta}({\bm{x}})))\}}{|B_{\delta}({\bm{x}})|},$ (2) $$

where $B_{\delta}({\bm{x}})$ is the ball centered at $x$ with a radius of $\delta$, and $C(N(\mathbb{R}\setminus B_{\delta}({\bm{x}})))$ represents the information about the point process outside the ball. If the Euclidean space is ordered, for instance, representing time, the conditioning term would represent the history of all points prior to $x$.

In general, effectively modeling and sampling from the conditional intensity (or related measures, e.g., hazard function or conditional density), for arbitrary metric spaces is generally not possible (Daley et al., ; Daley & Vere-Jones, ).
This difficulty has led to a variety of simplified parametrizations that restrict the captured point interactions (Ozaki, ; Zhou & Yu, ; Daley et al., ; Daley & Vere-Jones, ); discretizations of the space (Ogata, ; Osama et al., ); and numerical or Monte Carlo approximations (Chen et al., ; Hong & Shelton, ).

In contrast, we propose a method that bypasses the abstract concept of a (conditional) intensity function by directly manipulating point sets through a latent variable model. Our approach leverages the following point process properties:(^1^11We provide a proof of both properties for general Borel sets in A.2.)

Superposition: Given two point processes $N_{1}$ and $N_{2}$ with intensities $\lambda_{1}$ and $\lambda_{2}$ respectively, we define the superposition of the point processes as $N=N_{1}+N_{2}$ or equivalently $X_{1}\bigcup X_{2}$. Then, the resulting point process $N$ has intensity $\lambda=\lambda_{1}+\lambda_{2}$.
Independent thinning: Given a point process $N$ with intensity $\lambda$, randomly removing each point with probability $p$ is equivalent to sampling points from a point process with intensity $(1-p)\lambda$.

### 2.2 Diffusion Models

Ho et al. () and Sohl-Dickstein et al. () introduced a new class of generative latent variable models – probabilistic denoising diffusion models.
Conceptually, these models learn to reverse a probabilistic nosing process to generate new data and consist of three main components: a noising process, a denoising process, and a sampling procedure.
The noising process is defined as a forward Markov chain $q(X_{t+1}|X_{t})$, which progressively noises a data sample $X_{0}\sim p_{\mathrm{data}}(X)$ over $T$ steps, eventually transforming it into a sample from a stationary noise distribution $X_{T}\sim p_{\mathrm{noise}}(X)$.
Then, the denoising process is learned to reverse the noising process by approximating the posterior $q(X_{t}|X_{0},X_{t+1})$ with a model $p_{\theta}(X_{t}|X_{t+1})$.
Finally, the sampling procedure shows how to generate samples from the learned data distribution $p_{\theta}(X)=\int p_{\mathrm{noise}}(X_{T})\prod_{t=0}^{T-1}p_{\theta}(X_{t}|
X_{t+1})\mathop{}\!\mathrm{d}X_{1}\dots\mathop{}\!\mathrm{d}X_{T}$.

## 3 Point Set Diffusion

In this section, we derive a diffusion-based latent variable model for point sets on general metric spaces by systematically applying the thinning and superposition properties of random sets.
This approach allows direct manipulation of random point sets, avoiding the need for the abstract concept of an intensity function.
We begin by outlining the forward noising process in Section 3.1, which stochastically interpolates between point sets from the generating process and those from a noise distribution.
Subsequently, we demonstrate how to learn to reverse this noising process to generate new random point sets in Section 3.2.
Finally, in Section 3.3, we show how to sample from our unconditional model and generate conditional samples for general conditioning tasks on the metric space.

### 3.1 Forward Process

Let $X_{0}\sim p_{\mathrm{data}}(X)$ be an i.i.d. sample from the generating point process, and let $X_{T}\sim p_{\mathrm{noise}}(X)$ represent a sample from a noise point process.
We define the forward process as a stochastic interpolation between the point sets $X_{0}$ and $X_{T}$ over $T$ steps.
This process is modeled as a Markov chain $q(X_{t+1}|X_{t})$, where $X_{t}$ is the superposition of two random subsets: $X_{t}^{\mathrm{thin}}\subseteq X_{0}$ and $X_{t}^{{\epsilon}}\subseteq X_{T}$.
Specifically, $\forall t:X_{t}=X_{t}^{\mathrm{thin}}\bigcup X_{t}^{{\epsilon}}$, where $X_{t}^{\mathrm{thin}}$ and $X_{t}^{{\epsilon}}$ are independent samples from a thinning and a noise process, respectively.
We define the thinning and noise processes given two noise schedules $\{\alpha_{t}\in(0,1)\}_{t=1}^{T}$ and $\{\beta_{t}\in(0,1)\}_{t=1}^{T}$ as follows:

Figure: Figure 2: The forward process is a Markov Chain $q(X_{t+1}|X_{t})$, that stochastically interpolates a data sample $X_{0}$ with a noise point set $X_{T}$ over $T$ steps by applying a thinning and a noise process.
Refer to caption: x2.png

Thinning Process: This process progressively thins points in $X_{0}^{\mathrm{thin}}=X_{0}$, removing signal over time.
At every step $t+1$, each point ${\bm{x}}\in X_{t}^{\mathrm{thin}}$ is independently thinned with probability $1-\alpha_{t+1}$:

$$ $q({\bm{x}}\in X_{t+1}^{\mathrm{thin}}|{\bm{x}}\in X_{t}^{\mathrm{thin}})= \alpha_{t+1}.$ (3) $$

Consequently, the thinning defines $n$ independent Bernoulli chains, and the probability of any point ${\bm{x}}\in X_{0}$ remaining in $X_{t}^{\mathrm{thin}}$ is:

$$ $q({\bm{x}}\in X_{t}^{\mathrm{thin}}|{\bm{x}}\in X_{0})=\bar{\alpha}_{t},$ (4) $$

where $\bar{\alpha}_{t}=\prod_{i=1}^{n}\alpha_{i}$.
Equivalently, the intensity of the thinned points at step $t$ is given by $\lambda_{t}^{\mathrm{thin}}=\bar{\alpha}_{t}\lambda_{0}$ and the number of remaining points follows a Binomial distribution: $Pr[n_{t}|X_{0}^{\mathrm{thin}}]=\text{Binomial}(|X_{0}|,\bar{\alpha}_{t}),
\text{where }n_{t}=|X_{t}^{\mathrm{thin}}|$.

Noise Process: This process adds random points $X_{T}^{{\epsilon}}\sim p_{\mathrm{noise}}(X)$ sampled from a noise point process with intensity $\lambda^{{\epsilon}}$.
At step $t+1$, we express $X_{t+1}^{{\epsilon}}|X_{t}^{{\epsilon}}$ as:

$$ $X_{t+1}^{{\epsilon}}=X_{t}^{{\epsilon}}\cup X_{t+1}^{\Delta{\epsilon}},\quad \text{where }X_{t+1}^{\Delta{\epsilon}}\sim\beta_{t+1}\lambda^{{\epsilon}}.$ (5) $$

By the superposition property, the intensity of $X_{t}^{{\epsilon}}$ is $\lambda_{t}^{{\epsilon}}=\bar{\beta}_{t}\lambda^{{\epsilon}},\text{where }\bar
{\beta}_{t}=\sum_{i=1}^{t}\beta_{i}$ and $\bar{\beta}_{t}\in[0,1]$.
Alternatively, we can view the noise process as a reversed thinning process: we sample $X_{T}^{{\epsilon}}\sim p_{\mathrm{noise}}(X)$ and thin it by $1-\bar{\beta}_{t}$ to obtain $X_{t}^{{\epsilon}}$.
Given a noise sample $X_{T}^{{\epsilon}}$, we then find that:

$$ $q({\bm{x}}\in X_{t}^{{\epsilon}}|{\bm{x}}\in X_{T}^{{\epsilon}})=\bar{\beta}_{ t}.$ (6) $$

Notably, this process is independent of the random point set $X_{0}$, i.e., $\forall t:q(X_{t}^{{\epsilon}}|X_{0})=q(X_{t}^{{\epsilon}}).$

We present a visual depiction of the two forward processes in Figure 2.
Finally, given that $\forall t:X_{t}=X_{t}^{\mathrm{thin}}\bigcup X_{t}^{{\epsilon}}$ it follows that for $\lim_{t\to T}\bar{\alpha}_{t}=0$ and $\lim_{t\to T}\bar{\beta}_{t}=1$ the stationary distribution is $q(X_{T}|X_{0})=p_{\mathrm{noise}}(X)$, which can be seen by applying the superposition property and finding the intensity of $X_{t}|X_{0}$ to be $\bar{\alpha}_{t}\lambda_{0}+\bar{\beta}_{t}\lambda^{{\epsilon}}$.
To summarize, the forward process gradually removes points from the original point set $X_{0}\sim p_{\mathrm{data}}(X)$ while progressively adding points of a noise point set $X_{T}\sim p_{\mathrm{noise}}(X)$, stochastically interpolating between data and noise.

### 3.2 Reverse Process

To generate samples from our diffusion model, i.e., $X_{T}\to\cdots\to X_{0}$, we need to learn how to reverse the forward process by approximating the posterior $q(X_{t}|X_{0},X_{t+1})$ with a model $p_{\theta}(X_{t}|X_{t+1})$.
We will start by deriving the posterior $q(X_{t}|X_{0},X_{t+1})$ from the forward process $q(X_{t+1}|X_{t})$ and then show how to parameterize and train $p_{\theta}(X_{t}|X_{t+1})$ to approximate the posterior.

Since the forward process consists of two independent processes (thinning and noise) and noticing that $X_{t+1}^{\mathrm{thin}}=X_{0}\bigcap X_{t+1}$ and $X_{t+1}^{{\epsilon}}=X_{t+1}\setminus X_{0}$, the posterior can be derived in two parts:

Thinning posterior:
Since all points in $X_{t+1}^{\mathrm{thin}}$ have been retained from $t=0$, it follows that $X_{t+1}^{\mathrm{thin}}\subseteq X_{t}^{\mathrm{thin}}$.
Then for each point in ${\bm{x}}\in X_{0}\setminus X_{t+1}^{\mathrm{thin}}$, we derive then posterior using Bayes’ theorem, applying Equation 3, Equation 4 and the Markov property:

$$ $\displaystyle q({\bm{x}}\in X_{t}^{\mathrm{thin}}|{\bm{x}}\notin X_{t+1}^{ \mathrm{thin}},{\bm{x}}\in X_{0})$ $\displaystyle=\frac{q({\bm{x}}\notin X_{t+1}^{\mathrm{thin}}|{\bm{x}}\in X_{t} ^{\mathrm{thin}})q({\bm{x}}\in X_{t}^{\mathrm{thin}}|{\bm{x}}\in X_{0})}{q(x \notin X_{t+1}^{\mathrm{thin}}|{\bm{x}}\in X_{0})}$ (7) $\displaystyle=\frac{(1-\alpha_{t+1})\bar{\alpha}_{t}}{(1-\bar{\alpha}_{t+1})}= \frac{\bar{\alpha}_{t}-\bar{\alpha}_{t+1}}{1-\bar{\alpha}_{t+1}}.$ (8) $$

Thus, we can sample $X_{t}^{\mathrm{thin}}$ by superposition of $X_{t+1}^{\mathrm{thin}}$ and thinning $X_{0}\setminus X_{t+1}^{\mathrm{thin}}$ with Equation 7.

Figure: Figure 3: The posterior reverses the stochastic interpolation of $X_{0}\to X_{T}$ of the forward process by adding back thinned points from the thinning process and thinning point added in the noise process.
Refer to caption: x3.png

Noise posterior: Following the reverse thinning interpretation of the noise process, each point in $X_{t}^{{\epsilon}}$ must have been in both $X_{t+1}^{{\epsilon}}$ and $X_{T}^{{\epsilon}}$.
Hence, we derive the posterior for each point in $X_{t+1}^{{\epsilon}}$ to still be in $X_{t}^{{\epsilon}}$ by following Equation 6, along with the fact that $X_{t}^{{\epsilon}}$ is independent from $X_{0}$:

$$ $\displaystyle q({\bm{x}}\in X_{t}^{{\epsilon}}|{\bm{x}}\in X_{t+1}^{{\epsilon} },{\bm{x}}\notin X_{0})$ $\displaystyle=\frac{q({\bm{x}}\in X_{t+1}^{\mathrm{thin}}|{\bm{x}}\in X_{t}^{ \mathrm{thin}})q({\bm{x}}\in X_{t}^{\mathrm{thin}})}{q({\bm{x}}\in X_{t+1}^{ \mathrm{thin}})}$ (9) $\displaystyle=\frac{1\cdot\bar{\beta}_{t}}{(\bar{\beta}_{t+1})}=\frac{\bar{ \beta}_{t}}{\bar{\beta}_{t+1}}.$ (10) $$

Thus, we can sample $X_{t}^{{\epsilon}}$ by thinning $X_{t+1}^{{\epsilon}}$ with probability $(1-\frac{\bar{\beta}_{t}}{\bar{\beta}_{t+1}})$.

#### Parametrization:

Given $X_{0}$ and $X_{T}$, the derived posterior can reverse the noising process to generate $X_{0}$.
However, to generate a new approximate sample $X_{0}\sim p_{\mathrm{data}}(X)$, we need to be able to sample from the posterior $q(X_{t}|X_{0},X_{t+1})$ without knowing $X_{0}$.
For this reason we approximate the posterior with a model $p_{\theta}(X_{t}|X_{t+1})$, where we choose $p_{\theta}(X_{t}|X_{t+1})=\int q(X_{t}|\widetilde{X}_{0},X_{t+1})p_{\theta}(
\widetilde{X}_{0}|X_{t+1})\mathop{}\!\mathrm{d}\widetilde{X}_{0}$ and training a neural network $p_{\theta}(\widetilde{X}_{0}|X_{t})$ to approximate $X_{0}|X_{t+1}$ for each $t+1$.

To effectively train this model, we have to condition our model $p_{\theta}(\widetilde{X}_{0}|X_{t})$ on $X_{t}$.
We propose to embed the points ${\bm{x}}\in X_{t}$ permutation invariant with a Transformer encoder with full attention and apply a sinusoidal embedding to embed $n=|X_{t}|$ and $t$.
Then, to probilistically predict $X_{0}|X_{t}$, we make use of the following case distinction for $X_{t}^{\mathrm{thin}}=X_{0}\bigcap X_{t}$ and $X_{0}\setminus X_{t}$:

First, predicting the retained points in $X_{t}$, i.e., the intersection of $X_{0}$ and $X_{t}$, is a binary classification task for which we train a multi-layer-perceptron (MLP) $g_{\theta}({\bm{x}}\in X_{t}^{\mathrm{thin}}|X_{t},t)$ with binary cross entropy loss $\mathcal{L}_{\mathrm{BCE}}$.
Second, the thinned points in $X_{0}$, i.e., $X_{0}\setminus X_{t}$, is a point set $N$, which can be represented by its counting measure, as a mixture of $n$ Dirac measures:

$$ $N=\sum_{i=1}^{n}\delta_{X_{i}}.$ (11) $$

In A.3, we prove that any finite mixture of Dirac deltas, such as $N$, can be approximated by an $L^{2}$ function in $L^{2}(D,\mu)$ for any metric space $D$.
In Euclidean spaces, we approximate the Dirac measure with a mixture of multivariate Gaussian distributions with diagonal covariance matrices.
Note that the multivariate Gaussian density function is a standard approximation of the Dirac delta function and, as the determinant of a diagonal covariance matrix $\bm{\Sigma}\coloneqq\sigma\bm{I}$ approaches zero, the Gaussian increasingly resembles the Dirac delta (See Section A.3). We parameterize the number of points to sample $n_{\theta}$ and the components of the mixture — weights $\bm{w}_{\theta}$, mean $\bm{\mu}_{\theta}$ and diagonal covariance matrix $\bm{\Sigma}_{\theta}$ — with an MLP $f_{\theta}$ and train it with the negative log likelihood $\mathcal{L}_{\mathrm{NLL}}$.

Lastly, to ensure the expected number of points at any time $t$ throughout the diffusion process is constant, we use $\bar{\alpha}_{t}=1-\bar{\beta}_{t}$ and a noise process with a constant intensity such that $\int_{A}\lambda^{{\epsilon}}=E[N(A)]$ for the bounded Borel set $A$ that represents our domain.

### 3.3 Sampling Procedure

Unconditional sampling:
Starting from a sample $X_{T}$ of the noise distribution, we apply our Point Set Diffusion model to sample a new $X_{0}$ over $T$ steps. We start by sampling $X_{T}\sim\lambda_{{\epsilon}}$ and then for all $t\in(T,\dots,1)$ sample $\widetilde{X}_{0}\sim p_{\theta}(X_{0}|X_{t})$ to subsequently apply the denoising posterior $q(X_{t-1}|\widetilde{X}_{0},X_{t})$ and attain $X_{t-1}$.
Finally, at step $1$ we sample $\widetilde{X}_{0}\sim p_{\theta}(X_{0}|X_{1})$. We present the extended sampling algorithm in Algorithm 2.

Figure: Figure 4: Examples of conditioning masks for $\mathbb{R}_{\geq 0}$ and $\mathbb{R}^{2}$.
Refer to caption: x4.png

Conditional sampling: Let $C:D\rightarrow\{0,1\}$ be a conditioning mask on our metric space $D$, where we define the masking of a subset $X\subseteq D$ as $C(X)\coloneq\{{\bm{x}}\in X|C({\bm{x}})=1\}$ and its complement as $C^{\prime}(X)\coloneq\{{\bm{x}}\in X|C({\bm{x}})=0\}$.
Then, we can leverage our Point Set Diffusion model to conditionally generate random point sets outside the conditioning mask by applying Algorithm 1:
Algorithm 1 Conditional sampling
1:$X_{0}^{c}=C(X_{0})$
2:$X_{T}\sim\lambda_{{\epsilon}}$
3:for $t=T,\dots,1$ do
4:     $\widetilde{X}_{0}\sim p_{\theta}(X_{0}|X_{t})$
5:     $\widetilde{X}_{t-1}\sim q(X_{t-1}|\widetilde{X}_{0},X_{t})$
6:     $X_{t-1}^{c}\sim q(X_{t-1}|X_{0}^{c})$
7:     $X_{t-1}=C^{\prime}(\widetilde{X}_{t-1})\cup C(X_{t-1}^{c})$
8:end for
9:return $C^{\prime}(X_{0})$
Thus, following this sampling procedure, we can generate conditional samples for any conditioning mask $C$, where we represent some illustrative conditioning masks for bounded sets on $\mathbb{R}_{\geq 0}$ and $\mathbb{R}^{2}$ depicting temporal forecasting, history prediction and general imputation tasks in Figure 4.

## 4 Experiments

Although point processes are fundamentally generative models, the standard evaluation method relies on reporting the negative log-likelihood (NLL) on a hold-out test set, effectively reducing the evaluation to future single-event predictions for STPPs and TPPS.
However, this approach presents two key issues.
First, computing the NLL depends on the implementation and parameterization of the (conditional) intensity function and is intractable for many models, necessitating approximations using Monte Carlo methods, numerical integration, or the evidence lower bound (ELBO), complicating a fair comparisons between models.
Second, evaluating the likelihood of each point conditioned on ground-truth points does not necessarily reflect how well a model captures the actual data distribution or its ability to perform on complex conditional generation tasks (Shchur et al., ).
To overcome these limitations, we evaluate the generative capabilities of our proposed Point Set Diffusion model by benchmarking it on a range of unconditional and conditional generation tasks for both SPP and STPP.
The training of our model and the hyperparameters are presented in A.5, while all baselines are trained reproducing their reported NLL following their proposed hyperparameters and code.

### 4.1 Data

We follow Chen et al. () and evaluate our model on four benchmark datasets: three real-world datasets — Japan Earthquakes (U.S. Geological Survey, ), New Jersey COVID-19 Cases (The New York Times, ), and Citibike Pickups (Citi Bike, ) —and one synthetic dataset, Pinwheel, based on a multivariate Hawkes process (Soni, ).
The pre-processing and splits of the datasets are identical to Chen et al. ().

### 4.2 Metrics

To evaluate both unconditional and conditional tasks, we compute distances between point process distributions and individual point sets, assuming the space is normed, and all points are bounded, i.e., $\forall i,\bm{x}_{i}\in[-1,1]^{d}$. We use the following metrics in our evaluation:

Sequence Length (SL):
To compare the length distribution of point sets, we report the Wasserstein distance between the two categorical distributions.
For conditional tasks, we compare the length of the generated point set to the ground truth by reporting the Mean Absolute Error (MAE).

Counting Distance (CD):
Xiao et al. () introduced a Wasserstein distance for ordered TPPs based on Birkhoff’s theorem.
We generalize this counting distance to higher-dimensional ordered Euclidean spaces (e.g., STPPs) using the $L_{1}$ distance:

$$ $CD(X,Y)=\frac{1}{d}\sum_{i=1}^{k}{||{\bm{x}}_{i}-{\bm{y}}_{i}||_{1}}+\sum_{j=k +1}^{l}{||U-{\bm{y}}_{j}||_{1}},$ (12) $$

Figure: Figure 5: SPP conditioning task: top ground truth, middle Regularized Method and bottom Point Set Diffusion.
Refer to caption: x5.png

where $X=\{{\bm{x}}_{i}\}_{i=1}^{k}$ and $Y=\{{\bm{y}}_{i}\}_{i=1}^{l}$ are two ordered samples from a point process on a metric space of dimensionality $d$, i.e. $D\subseteq\mathbb{R}^{d}$. Further, $U\coloneqq({\bm{u}}_{1},\dots,{\bm{u}}_{d})$ represents the upper bounds of the metric space $D$ along each dimension and we assume, without loss of generality, $l\geq k$.

Wasserstein Distance (WD):
An instance of a Point Process is itself a stochastic process of points in space. Hence, we can compute a distance between two point sets based on the Wasserstein distance on the metric space $D\subseteq\mathbb{R}^{d}$ between the two sets of points.

Maximum Mean Discrepancy (MMD) (Gretton et al., ):
The kernel-based statistic test compares two distributions based on a distance metric; we use the WD for SPPs and CD for STPP.

**Table 1: Density estimation results on the hold-out test set for SPPs, averaged over three random seeds (bold best and underline second best).**
|  | Earthquakes | Covid NJ | Citybike | Pinwheel |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | SL($\downarrow)$ | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) |
| Log-Gaussian Cox | $\underline{0.047}$ | $\underline{0.214}$ | $\underline{0.209}$ | $\underline{0.340}$ | $0.104$ | $\underline{0.336}$ | 0.017 | $\underline{0.285}$ |
| Regularized Method | $2.361$ | $0.391$ | $0.255$ | $0.411$ | $\underline{0.097}$ | $0.342$ | $\underline{0.039}$ | $0.411$ |
| Point Set Diffusion | 0.038 | 0.173 | 0.199 | 0.268 | 0.056 | 0.092 | 0.017 | 0.099 |

### 4.3 Spatial Point Processes

We evaluate our model’s ability to capture the distribution of spatial point processes (SPP) by benchmarking it against two methods. The first is the widely used Log-Gaussian Cox Process (Jesper Møller, ), a doubly stochastic model that parameterizes the intensity function using a Gaussian process.
The second is the Regularized Method (Osama et al., ), a spatial model that leverages a regularized criterion to infer predictive intensity intervals, offering out-of-sample prediction guarantees and enabling conditional generation.

Unconditional Generation (Density Estimation):
In this experiment, we generate 1,000 unconditional samples from each model and compare their distribution to a hold-out test set using the WD-SL and WD-MMD metrics.
As shown in Table 2, our Point Set Diffusion model consistently generates samples most closely matching the data distribution across all datasets.
While the baseline models perform reasonably well in capturing the count distributions for most datasets, their reliance on spatial discretization and smoothness properties of the intensity function limit their ability to capture the complex spatial patterns in the data, as reflected by higher WD-MMD scores.

Conditional Generation: To assess Point Set Diffusion’s ability to solve spatial conditioning tasks, we sample 50 random bounding boxes (with widths uniformly sampled between 1/8 and 3/8 of the metric space) for imputation on the hold-out test set, and report the results in Table 2.
The Regularized Method fits a spatial Poisson model with out-of-sample accuracy guarantees and has been shown by Osama et al. () to outperform the Log-Gaussian Cox Process on interpolation and extrapolation tasks.
However, we find that the Regularized Method’s reliance on predicting a smooth and discretized intensity function conditioned on neighboring areas leads to inaccurate imputations when the adjacent regions contain significantly different numbers of points (see hexagonal discretization structure and smoothness in Figure 5).
This issue is exacerbated by not capturing a shared intensity function across point sets, making it difficult for the Regularized Method to handle non-smooth spatial patterns, such as varying inhomogeneous intensities shared across multiple point sets.
This highlights a core limitation of SPP models that rely on instance-specific intensity functions.

### 4.4 Spatio-temporal Point Processes

For STPPs, we evaluate our model’s ability to capture the point process distribution by benchmarking it against three state-of-the-art STPP models, learning an autoregressive intensity function.

Figure: Figure 6: STPP runtime for sampling $n$ points.
Refer to caption: x6.png

DeepSTPP (Zhou et al., ) uses a latent variable framework to non-parametrically model the conditional intensity based on kernels.
DiffSTPP (Yuan et al., ) is based on a diffusion model approximating the conditional intensity.
Lastly, AutoSTPP (Zhou & Yu, ) uses the automatic integration for neural point processes, presented by Lindell et al. (), to parameterize a generalized spatiotemporal Hawkes model.

Sampling Runtime: We report the median sampling runtime over ten runs generating ten point set of length $n$ on an NVIDIA A100-PCIE-40GB for all STPP models in Figure 6. Point Set Diffusion achieves a near constant sampling runtime for all point set lengths as it generates all points in parallel, while all autoregressive baselines, given their sequential sampling, show at least a linear relationship between runtime and $n$.

Unconditional Generation (Density Estimation):
We evaluate the performance of each model by comparing the WD-SL and CS-MMD between the hold-out test set and 1,000 samples generated by the trained models, as shown in Table 4.
Again, the Point Set Diffusion model best captures the distribution of the point process distribution for all datasets.
The autoregressive intensity functions of the baseline models fail to generate point sets that align closely with the data distribution for most datasets, as reflected in the stark differences in the WD-SL and CD-MMD metrics compared to Point Set Diffusion.
While these baselines are trained to predict the next event given a history window, they struggle to unconditionally sample realistic point sets when starting from an empty sequence.
Consequently, this highlights our argument that the standard evaluation based on NLL is insufficient to assess the true generative capacity of point process models.

Conditional Generation (Forecasting):
Forecasting future events based on historical data is a challenging and a fundamental task for STPP models.
To evaluate this capability, we uniformly sampled 50 random starting times from the interval $[\frac{5}{8}U_{\text{time}},\frac{7}{8}U_{\text{time}}]$, where $U_{\text{time}}$ is the maximum time, for each point set in the hold-out test set.
The results are detailed in Table 4.(^2^22AutoSTPP is not included in this analysis due to its prohibitively slow sampling speed (see Figure 6 and the limitations discussed in Zhou & Yu ()), which made it impractical to sample the 50 forecast windows for all instances in the test set within a reasonable timeframe.)
The autoregressive baselines, trained to predict the next event based on history, achieve good forecasting results for most datasets, one even surpassing Point Set Diffusion on the Covid NJ dataset.
Still, our unconditional model outperforms the autoregressive baselines across all other datasets.

**Table 3: Density estimation results on the hold-out test set for STPP, averaged over three random seeds (bold best and underline second best).**
|  | Earthquakes | Covid NJ | Citybike | Pinwheel |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) |
| DeepSTPP | 1.489 | 0.232 | 0.300 | 0.439 | 2.403 | 0.327 | 0.470 | 0.167 |
| DiffSTPP | 0.088 | 0.064 | 0.332 | 0.146 | 0.560 | 0.611 | 0.196 | 0.055 |
| AutoSTPP | 0.169 | 0.472 | 0.390 | 0.382 | 0.556 | 0.412 | 0.111 | 0.105 |
| Point Set Diffusion | 0.042 | 0.023 | 0.189 | 0.043 | 0.032 | 0.020 | 0.023 | 0.020 |

### 4.5 Other conditioning task

Since the STPP baselines are auto-regressive models, they are limited to forecasting tasks.
However, our model can generate conditional samples for any conditioning mask $C$ on our metric space.
To showcase this feature of our model, we present a few visual examples of complex conditioning tasks in Figure 7.

## 5 Related Work

Since large parts of the real-world can be effectively captured by Euclidean spaces, point processes have mainly been defined on spatial and temporal dimensions, represented by an Euclidean space.
Hence, for this discussion of the related work, we will focus on unordered and ordered point processes on Euclidean spaces, mainly SPPs, TPPs and STPPs.
For completeness, we want to mention traditional parametric point processes defined on manifolds, such as determential point processes (Berman, ; Katori & Shirai, ) and cluster point processes (Bogachev & Daletskii, ).

Unordered Point Processes (SPP):
Modeling a permutation-invariant intensity for unordered point sets that captures complex interactions while remaining efficient for sampling is challenging (Daley & Vere-Jones, ), seemingly limiting the development of machine-learning-based models for SPPs.
Classical models like the Poisson Point Process (Kingman, ) use either homogeneous or inhomogeneous intensity functions across space.
More flexible models, such as Cox processes (Cox, ), and specifically the popular Log-Gaussian Cox Process (Jesper Møller, ), extend this by modeling the intensity function through a doubly stochastic process, allowing for flexible spatial inhomogeneity.
A recent approach, the Regularized Method by Osama et al. (), parameterizes a spatial Poisson process on a hexagonal grid with splines, offering out-of-sample guarantees.
However, these methods often rely on spatial discretization and simple parametric forms and some require separate intensity estimates for each point set, limiting their ability to capture the underlying distribution across different samples (Daley & Vere-Jones, ; Osama et al., ).

Figure: Figure 7: Complex spatial conditioning tasks solved with Point Set Diffusion: Top condition and ground truth data, bottom density plots for predictions.
Refer to caption: x7.png

Ordered point processes (TPP and STPP):
The causal ordering of time enables the parametrization of a conditional intensity, which classically is being modeled with parametric functions, where the Hawkes Process (Hawkes, ) is the most widely used model and captures point interaction patterns like self-excitation.
Given the sequential nature of ordered point process a variety of Machine Learning based approaches for TPPs and STPPs have been proposed (see Shchur et al. () for a review on neural TPPs).
Where, recurrent neural network- (Du et al., ; Shchur et al., ) and transformer-based encoders (Zhang et al., ; Zuo et al., ; Chen et al., ) are leveraged to encode the history and neurally parameterized Hawkes (Zhou & Yu, ; Zhang et al., ; Zuo et al., ), parametric density functions (Du et al., ; Shchur et al., ), mixtures of kernels (Okawa et al., ; Soen et al., ; Zhang et al., ; Zhou et al., ), neural networks (Omi et al., ; Zhou & Yu, ), Gaussian diffusion (Lin et al., ; Yuan et al., ) and normalizing flows (Chen et al., ; Shchur et al., ) have been proposed to (non)-parametrically decode the conditional density or intensity of the next event.

Differences to Add-Thin (Lüdke et al., ):
Since our method is closely related to Add-Thin, we want to highlight their key methodological differences.
While Add-Thin proposed to leverage the thinning and superposition properties to define a diffusion process for TPPs, Point Set Diffusion generalizes this idea to define a diffusion-based latent variable model for point processes on general metric spaces.
In doing so, we disentangle the superposition and thinning to attain two independent processes to allow for more explicit control and define the diffusion model independent of the intensity function as a stochastic interpolation of point sets.
Furthermore, Add-Thin has to be trained for specific conditioning tasks, while we show how to condition our unconditional Point Set Diffusion model for arbitrary conditioning tasks on the metric space.
Lastly, Point Set Diffusion and its parametrization is agnostic to the ordering of points, making it applicable to model the general class of point processes on any metric space, including for example SPPs.

## 6 Conclusion

To model general point processes on metric spaces, we present Point Set Diffusion, a novel diffusion-based latent variable model.
We derive Point Set Diffusion as a stochastic interpolation between data point sets and noise point sets governed by the thinning and superposition properties of random point sets.
Thereby, we attain a very flexible, unconditional Point Process model that can be conditioned for arbitrary condition masks on the metric space and allows for efficient and parallel sampling of entire point sets without relying on the (conditional) intensity function.
In conditional and unconditional experiments on synthetic and real-world SPP and STPP data, we demonstrate that Point Set Diffusion achieves state-of-the-art performance while allowing for up to orders of magnitude faster sampling.

With Point Set Diffusion, we have presented a novel set modeling approach and would be interested to see how future work explores it’s limitations on other (high-dimensional) metric (e.g., Riemannian manifolds), topological and discrete spaces with potential applications extending beyond traditional point sets including but not limited to natural language and graphs.

## Appendix A Appendix

### A.1 Sampling Algorithm

Figure: Algorithm 2 Sampling

### A.2 Point Process Properties

The thinning and superposition properties have been proved by other works for different versions of point processes. For completeness and generality, we prove them for a general Borel set $A$. To apply these proofs for SPPs consider $A\subseteq\mathcal{S}$, where $\mathcal{S}$ is a metric space in $\mathbb{R}^{d}$ and for STPPs consider $A\subseteq[0,T]\times\mathcal{S}$, where $T>0$.

#### Superposition:

Proof.
It is straightforward to obtain the superposition expectation measure from Equation 1:

$$ $\mu(A)=\mathbb{E}[N(A)]=\mathbb{E}[N_{1}(A)+N_{2}(A)]=\mathbb{E}[N_{1}(A)]+ \mathbb{E}[N_{2}(A)]=\mu_{1}(A)+\mu_{2}(A).$ (13) $$

Then, every point process has an intensity of $\lambda_{1}$ and $\lambda_{2}$ for each of the expectation measures $\mu_{1}$ and $\mu_{2}$, respectively. Therefore, taking the right-hand side of Equation 1, we obtain the following intensity function for the superposition of point processes:

$$ $\mu(A)=\mu_{1}(A)+\mu_{2}(A)=\int_{A}{\lambda_{1}(x)dx}+\int_{A}{\lambda_{2}(x )dx}=\int_{A}{\lambda_{1}(x)+\lambda_{2}(x)dx}.$ (14) $$

This states that the density function for expectation measure $\mu$ is $\lambda\coloneqq\lambda_{1}+\lambda_{2}$, and concludes the proof for the superposition property of intensities for point processes.
$\square$

#### Thinning:

Proof.
For this property, we need to assume that the singletons are simple, so we can only have one point at each position: $N(\{x\})\leq 1$; these point processes are called simple. Simple point processes can be represented as a sum of Dirac measures at the random points $X_{i}\in\mathcal{S}$:

$$ $N=\sum_{i}{\delta_{X_{i}}}.$ (15) $$

The previous assumption on singletons makes the sum above a finite sum. If $Z_{i}\in\{0,1\}$ are Bernoulli random variables with a success probability $p$ we can define a random thinning process as the superposition of the following point processes:

$$ $N_{1}=\sum_{i}{Z_{i}\delta_{X_{i}}}.$ (16) $$

$$ $N_{2}=\sum_{i}{(1-Z_{i})\delta_{X_{i}}}.$ (17) $$

Since there are only two options for the Bernoulli random variable, it holds that the superposition of the point processes defined in Equations 16 and 17 are equivalent to the original point process, i.e., $N=N_{1}+N_{2}$.
Given that all $Z_{i}\sim Bern(p)$ are i.i.d., we obtain a conditional probability distribution on the thinned point process $N_{1}(A)|N(A)=n\sim Binom(n,p).$ And by the law of total expectation, we derive:

$$ $\mu_{1}(A)=\mathbb{E}[N_{1}(A)]=\mathbb{E}\big{[}\mathbb{E}[N_{1}(A)|N(A)]\big {]}=\mathbb{E}[N(A)p]=\mu(A)p.$ (18) $$

We can write the terms of the equation before in terms of the intensity measure of the point process:

$$ $\mu_{1}(A)=p\cdot\mu(A)=p\int_{A}{\lambda(x)dx}=\int_{A}{p\lambda(x)dx}.$ (19) $$

Hence, the equation above implies that the intensity of the new point process $N_{1}$, which keeps the points of the original point process $N$ with probability $p$, is $p\lambda$.
By the property of superposition, since $N=N_{1}+N_{2}$, then $\lambda=p\lambda+(1-p)\lambda$. Therefore, the intensity of the point process $N_{2}$,
containing the thinned points, is $(1-p)\lambda$.
This proves that, in the opposite case, when removing points with probability $p$ from a given point process with intensity $\lambda$, the intensity of the point process with the points kept after thinning is $(1-p)\lambda$.
$\square$

### A.3 Approximation of Mixture of Dirac Delta Functions by L 2 superscript 𝐿 2 L^{2} italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -Functions

###### Definition 1 (Dirac delta function)

Let $(D,d,\mu)$ be a general metric space equipped with a measure $\mu$. A Dirac delta function $\delta_{\bm{x}}$ at a point $\bm{x}\in D$ is defined as a distribution such that for any test function $f$:

$$ $\int_{D}{f(\bm{y})\delta_{\bm{x}}(\bm{y})d\mu(\bm{y})}=f(\bm{x}).$ (20) $$

###### Theorem 1

Let $f_{M}(\bm{y})$ be a finite mixture of Dirac deltas:

$$ $f_{M}(\bm{y})=\sum_{i=1}^{n}w_{i}\delta_{\bm{x}_{i}}(\bm{y}),$ (21) $$

where $\bm{x}_{1},\dots,\bm{x}_{n}\in D$ are points in the metric space, and $w_{i}\in\mathbb{R}$ are weights associated with each Dirac delta function. Then, this finite mixture of Dirac deltas $f_{M}$ can be approximated by $L^{2}$ functions in $L^{2}(D,\mu)$.

Proof.
We use a sequence of smooth functions that approximate each Dirac delta in the mixture and then show that this approximation converges in the $L^{2}$-norm.

Firstly, we show how to approximate Dirac delta functions. Let us consider a family of smooth functions $\phi_{\epsilon}(\bm{x})$ (such as bump functions or mollifiers) that approximate the Dirac delta function $\delta_{\bm{x}}$ as $\epsilon\to 0$. These functions $\phi_{\epsilon}(\bm{x}-\bm{x}_{i})$ are supported near $\bm{x}_{i}$ and satisfy:

$$ $\lim_{\epsilon\to 0}\phi_{\epsilon}(\bm{x}-\bm{x}_{i})=\delta_{\bm{x}_{i}}(\bm {x}).$ (22) $$

In particular, for any test function $f$, we have:

$$ $\int_{D}{f(\bm{y})\phi_{\epsilon}(\bm{y}-\bm{x}_{i})d\mu(\bm{y})}\to f(\bm{x}_ {i})\quad\text{as}\quad\epsilon\to 0.$ (23) $$

Hence, $\phi_{\epsilon}(\bm{x}-\bm{x}_{i})$ has a similar property as the one of Dirac deltas given in Equation 20 and serves as an approximation of the Dirac delta $\delta_{\bm{x}_{i}}(\bm{x})$ for a small $\epsilon$.

Secondly, we approximate the mixture of Dirac deltas $f_{M}$ by a function in $L^{2}(D,\mu)$ using the same $\phi_{\epsilon}(\bm{x})$-based approximation for each Dirac delta, defining:

$$ $f_{\epsilon}(\bm{y})=\sum_{i=1}^{n}w_{i}\phi_{\epsilon}(\bm{y}-\bm{x}_{i}).$ (24) $$

Each term $\phi_{\epsilon}(\bm{y}-\bm{x}_{i})$ is a smooth approximation of the corresponding Dirac delta $\delta_{\bm{x}_{i}}(\bm{y})$, and the sum represents the approximation of the entire mixture of Diracs.

Thirdly, we show that the sequence $f_{\epsilon}$ converges to $f_{M}$ in the $L^{2}$-norm, i.e., that:

$$ $\lim_{\epsilon\to 0}\|f_{\epsilon}-f_{M}\|_{L^{2}(D,\mu)}=0.$ (25) $$

Since $f_{M}$ is a sum of Dirac deltas, it is not directly in $L^{2}(D,\mu)$, but its approximation $f_{\epsilon}$ is because each $\phi_{\epsilon}$ is a smooth function and smooth functions with compact support are in $L^{2}(D,\mu)$.

We compute now the squared $L^{2}$ norm of the difference in Equation 25:

$$ $\|f_{\epsilon}-f_{M}\|_{L^{2}(D,\mu)}^{2}=\int_{D}|f_{\epsilon}(\bm{y})-f_{M}( \bm{y})|^{2}\,d\mu(\bm{y}).$ (26) $$

Note that the squared difference of $f_{\epsilon}$ and $f_{M}$ in the above equation will have quadratic and crossed terms. However, we can neglect the crossed terms: $2\sum_{i<j}w_{i}w_{j}\int_{D}\left(\phi_{\epsilon}(\bm{y}-\bm{x}_{i})-\delta_{
\bm{x}_{i}}(\bm{y})\right)\left(\phi_{\epsilon}(\bm{y}-\bm{x}_{j})-\delta_{\bm
{x}_{j}}(\bm{y})\right)\,d\mu(\bm{y})$, since every smooth function $\phi_{\epsilon}(\bm{y}-\bm{x}_{i})$ is concentrated near $\bm{x}_{i}$ and terms involving different indices do not contribute to the limit.

Hence, we can simplify the norm in Equation 26 into the sum of the individual terms:

$$ $\|f_{\epsilon}-f_{M}\|_{L^{2}(D,\mu)}^{2}=\sum_{i=1}^{n}\int_{D}{w_{i}^{2}|( \phi_{\epsilon}(\bm{y}-\bm{x}_{i})-\delta_{\bm{x}_{i}}(\bm{y}))|^{2}d\mu(\bm{y })}.$ (27) $$

For every $i$, the term $\int_{D}{|\phi_{\epsilon}(\bm{y}-\bm{x}_{i})-\delta_{\bm{x}_{i}}(\bm{y})|^{2}d
\mu(\bm{y})}$ becomes small as $\epsilon\to 0$, because by construction $\phi_{\epsilon}(\bm{y}-\bm{x}_{i})\to\delta_{\bm{x}_{i}}(\bm{y})$ in the sense of distributions.
Thus, by the properties of $\phi_{\epsilon}$, we conclude that:

$$ $\lim_{\epsilon\to 0}\|f_{\epsilon}-f_{M}\|_{L^{2}(D,\mu)}=0.$ (28) $$

$\square$

###### Lemma 1

Let $p(\bm{x};\bm{\mu},\bm{\Sigma})$ be the probability density function (PDF) of a multivariate Gaussian distribution. Then $p\in L^{2}(\mathbb{R}^{d})$.

Proof.
The PDF of a multivariate Gaussian distribution in $\mathbb{R}^{d}$ with mean vector $\bm{\mu}\in\mathbb{R}^{d}$ and covariance matrix $\bm{\Sigma}$ (which is positive definite) is given by:

$$ $p(\bm{x};\bm{\mu},\bm{\Sigma})=\frac{1}{(2\pi)^{n/2}|\bm{\Sigma}|^{1/2}}\exp \left(-\frac{1}{2}(\bm{x}-\bm{\mu})^{T}\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu})\right),$ (29) $$

where $\bm{x}\in\mathbb{R}^{d}$, $|\bm{\Sigma}|$ is the determinant of the covariance matrix $\bm{\Sigma}$, and $\bm{\Sigma}^{-1}$ is the inverse of the covariance matrix. We show that $\|p\|_{L^{2}}=\left(\int_{\mathbb{R}^{d}}|p(\bm{x})|^{2}\,d\bm{x}\right)^{1/2}$ is finite.

We need to compute the following integral:

$$ $\int_{\mathbb{R}^{d}}p(\bm{x})^{2}\,d\bm{x}=\frac{1}{(2\pi)^{n}|\bm{\Sigma}|} \int_{\mathbb{R}^{d}}\exp\left(-(\bm{x}-\bm{\mu})^{T}\bm{\Sigma}^{-1}(\bm{x}- \bm{\mu})\right)d\bm{x}.$ (30) $$

To simplify the calculation, we perform a change of variables: $\bm{y}=\bm{\Sigma}^{-1/2}(\bm{x}-\bm{\mu})$. Under this transformation: $(\bm{x}-\bm{\mu})^{T}\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu})=\bm{y}^{T}\bm{y}=\|\bm{
y}\|^{2}$, and the differential $d\bm{x}$ transforms as: $d\bm{x}=|\bm{\Sigma}^{1/2}|\,d\bm{y}=|\bm{\Sigma}|^{1/2}\,d\bm{y}$. Substituting these into the integral, we get:

$$ $\int_{\mathbb{R}^{d}}\exp\left(-(\bm{x}-\bm{\mu})^{T}\bm{\Sigma}^{-1}(\bm{x}- \bm{\mu})\right)d\bm{x}=|\bm{\Sigma}|^{1/2}\int_{\mathbb{R}^{d}}\exp(-\|\bm{y} \|^{2})\,d\bm{y}=\pi^{n/2},$ (31) $$

since the remaining integral is a standard Gaussian integral. Thus, the $L^{2}$-norm integral becomes:

$$ $\int_{\mathbb{R}^{d}}p(\bm{x})^{2}\,d\bm{x}=\frac{1}{(2\pi)^{n}|\bm{\Sigma}|}| \bm{\Sigma}|^{1/2}\pi^{n/2}=\frac{1}{2^{n}\pi^{n/2}|\bm{\Sigma}|^{1/2}}.$ (32) $$

Since the integral is a finite constant, we conclude that the PDF belongs to $L^{2}(\mathbb{R}^{d})$.

$\square$

###### Corollary 1

Given an Euclidean space $D\subseteq\mathbb{R}^{d}$, a finite sum of Dirac deltas can be approximated with a mixture of multivariate Gaussian distributions:

$$ $p_{M}(\bm{x})=\sum_{i=1}^{n}{w_{i}\cdot\mathcal{N}(\bm{x};\bm{\bm{\mu}_{i}}, \bm{\Sigma}_{i})}.$ (33) $$

Proof.
Note that we do not show that a mixture of multivariate Gaussian distributions is the best candidate to approximate a finite sum of Dirac deltas. However, note that a multivariate Gaussian distribution is a standard approximation of a Dirac delta function and can, in the limit of a small covariance matrix, i.e. $|\bm{\Sigma}|<<1$, approximate it.

The aim of this proof is to show that the mixture of Gaussians $p_{M}$ can be a candidate to approximate the Dirac deltas. From Theorem 1, this is equivalent to showing that $p_{M}$ is a $L^{2}$ function.

To prove this, we need to integrate:

$$ $\int_{D}{p_{M}(\bm{x})^{2}d\bm{x}}=\sum_{i=1}^{k}w_{i}^{2}\int_{D}{\mathcal{N} (\bm{x};\bm{\mu}_{i},\bm{\Sigma}_{i})^{2}d\bm{x}}+2\sum_{i\neq j}w_{i}w_{j} \int_{D}{\mathcal{N}(\bm{x};\bm{\mu}_{i},\bm{\Sigma}_{i})\mathcal{N}(\bm{x}; \bm{\mu}_{j},\bm{\Sigma}_{j})d\bm{x}}.$ (34) $$

On the one hand, Lemma 1 shows that the integrals of the first sum are finite constants. On the other hand, the integrals on the second sum cannot be computed in a closed form, but it is well known that the product decays exponentially as $\|\bm{x}\|\to\infty$ ensuring a finite integral. Therefore, the squared integral is just a sum of finite constants and hence finite.
$\square$

###### Definition 1 (Dirac delta function)

###### Theorem 1

###### Lemma 1

###### Corollary 1

### A.4 Metrics with standard deviations

**Table 5: Density estimation results on the hold-out test set for SPPs averaged over three random seeds.**
|  | Earthquakes | Covid NJ | Citybike | Pinwheel |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | SL($\downarrow)$ | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) |
| Log-Gaussian Cox | 0.047$\pm$ 0.014 | 0.214$\pm$ 0.004 | 0.209$\pm$ 0.011 | 0.340$\pm$ 0.008 | 0.104$\pm$ 0.017 | 0.336$\pm$ 0.014 | 0.017$\pm$ 0.004 | 0.285$\pm$ 0.004 |
| Regularized Method | 2.361$\pm$ 0.064 | 0.391$\pm$ 0.004 | 0.255$\pm$ 0.011 | 0.411$\pm$ 0.003 | 0.097$\pm$ 0.008 | 0.342$\pm$ 0.008 | 0.039$\pm$ 0.003 | 0.411$\pm$ 0.004 |
| Point Set Diffusion | 0.038$\pm$ 0.003 | 0.173$\pm$ 0.004 | 0.199$\pm$ 0.002 | 0.268$\pm$ 0.016 | 0.056$\pm$ 0.020 | 0.092$\pm$ 0.020 | 0.017$\pm$ 0.003 | 0.099$\pm$ 0.006 |

**Table 7: Density estimation results on the hold-out test set for STPP averaged over three random seeds.**
|  | Earthquakes | Covid NJ | Citybike | Pinwheel |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) | SL($\downarrow$) | MMD($\downarrow$) |
| DeepSTPP | 1.489$\pm$ 0.988 | 0.232$\pm$ 0.046 | 0.300$\pm$ 0.093 | 0.439$\pm$ 0.257 | 2.403$\pm$ 0.884 | 0.327$\pm$ 0.040 | 0.470$\pm$ 0.326 | 0.167$\pm$ 0.032 |
| DiffSTPP | 0.088$\pm$ 0.009 | 0.064$\pm$ 0.024 | 0.332$\pm$ 0.012 | 0.146$\pm$ 0.026 | 0.560$\pm$ 0.045 | 0.611$\pm$ 0.113 | 0.196$\pm$ 0.098 | 0.055$\pm$ 0.005 |
| AutoSTPP | 0.169$\pm$ 0.047 | 0.472$\pm$ 0.275 | 0.390$\pm$ 0.020 | 0.382$\pm$ 0.107 | 0.556$\pm$ 0.077 | 0.412$\pm$ 0.159 | 0.111$\pm$ 0.019 | 0.105$\pm$ 0.005 |
| Point Set Diffusion | 0.042$\pm$ 0.003 | 0.023$\pm$ 0.003 | 0.189$\pm$ 0.006 | 0.043$\pm$ 0.003 | 0.032$\pm$ 0.004 | 0.020$\pm$ 0.001 | 0.023$\pm$ 0.003 | 0.020$\pm$ 0.001 |

### A.5 Model Setup

#### Architecture:

The classifier to predict $X_{0}\cap X_{t}$ is a MLP with $3$ layers and ReLU as activation function. The mixture of multivariate Gaussian distribution that approximates $X_{0}\setminus X_{t}$ contains $16$ components, and the parameters are learned with an MLP of $2$ layers and ReLU as an activation function.

#### Training:

All models have been trained on an NVIDIA A100-PCIE-40GB. We use Adam as the optimizer and a fixed weight decay of $0.0001$ to avoid overfitting. To avoid exploding gradients, we clip the gradients to have a norm lower than $2$.

#### Hyperparameters:

We use the same hyperparameters for all datasets and types of point processes. We leverage a hidden dimension and embedding size of $32$. For training, we use a batch size is chosen of $128$ and a learning rate of $0.001$.

#### Early stopping:

We train the models up to $5000$ epochs with early stopping, sampling $100$ sequences from the model and comparing them to the validation split, with WD-SL metric for SPP and the CD-MMD metric for STPPs.

### A.6 Connecting the applied loss and the ELBO

In this section we connect our applied loss function and the ELBO to the unknown data distribution. The ELBO of diffusion models Ho et al. () is given by:

$$ $\begin{split}\mathcal{L}_{ELBO}&=\mathbb{E}_{q}\Big{[}D_{KL}\big{(}q(X_{T}|X_{ 0})\,||\,p(X_{T})\big{)}-\log{p_{\theta}(X_{0}|X_{1})}\\ &+\sum_{t=2}^{T}{D_{KL}\Big{(}q(X_{t-1}|X_{0},X_{t})\,||\,p_{\theta}(X_{t-1}|X _{t})\Big{)}}\Big{]}.\end{split}$ (35) $$

The first term in the ELBO is constant as the distributions defining $q(X_{T}|X_{0})$ and $p_{noise}(X_{T})$ have not any learning parameters. The second term is minimized as we directly train our model to optimize this likelihood. The last term, consisting on a sum of different KL divergences between two densities, can be defined as:

$$ $D_{KL}(q\,||\,p_{\theta})=\mathbb{E}_{q}\big{[}\log{(q(X_{t-1}|X_{0},X_{t})}- \log{(p_{\theta}(X_{t-1}|X_{t}))}\big{]}.$ (36) $$

Hence, we minimize the KL divergence by minimizing the only $\theta$-dependent term, on the right-hand side, maximizing the expectation over the log-likelihood $\log{(p_{\theta}(X_{t-1}|X_{t}))}$. As mentioned in Section 3.2, we can divide the reverse process in two parts: the thinning posterior and the noise posterior.

#### Thinning posterior.

The component $X_{t-1}^{thin}\coloneqq X_{0}\cap X_{t-1}$ is composed of points kept from $X_{t}^{thin}$ and points thinned from $X_{0}\setminus X_{t}^{thin}$.

For retained points, both the posteriors $q(X_{t-1}|X_{0},X_{t})$ and $p_{\theta}(X_{t-1}|X_{t})$ are defined by Bernoulli distributions. The cross-entropy $H(q,p)$ can be written by definition as $H(q,p)=H(q)+D_{KL}(q\,||\,p)$, where $H(q)$ is the entropy. Therefore, minimizing the KL divergence is equivalent to minimizing the (binary) cross entropy.

For the other points, we sample from the posterior $q(X_{t-1}|X_{0},X_{t})$ by thinning $X_{0}\setminus X_{t}^{thin}$. For this, we approximate first $X_{0}$ as $X_{0}=X_{t}\cup(X_{0}\setminus X_{t})$, where the component $X_{0}\setminus X_{t}$ is learned by maximizing the log-likelihood $\mathbb{E}_{q}[\log{(p_{\theta}(X_{t-1}|X_{t}))}].$ Optimizing this term is equivalent to minimizing the negative log-likelihood (NLL) of our model’s mixture of multivariate Gaussian, representing the probability density of $X_{0}\setminus X_{t}$ through a mixture model.

#### Noise posterior.

The posterior $X_{t-1}^{\epsilon}$ can be derived by thinning $X_{t}^{\epsilon}=X_{t}\setminus X_{t}^{thin}$. In this case the posteriors $q(X_{t-1}|X_{0},X_{t})$ and $p_{\theta}(X_{t-1}|X_{t})$ directly depend on $X_{t}^{thin}$, which we approximate with the thinning’s posterior explained above by minimizing the cross entropy and, equivalently, minimizing the KL divergence.

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 5
presentation: 4
questions: 1. In Figure 5, can the authors show the predicted density of different trajectories in the same masked area?
soundness: 3
strengths: The idea of using the diffusion-style model to characterize point processes is super interesting. The content is clear and well-written, making the methodology and the results accessible to the reader. The paper also covers unconditional and conditional sampling methods, which have the potential to correspond to two important questions in the point process modeling (first-order and second-order modeling). The authors also provide thorough experimentation to validate the effectiveness of the proposed model.
summary: This paper proposes a novel modeling approach to point processes via diffusion on point sets (discrete events), addressing the reliance on the intensity function when establishing or learning the model. It can capture the distribution of point processes and generate a series of events based on noise point sets. Meanwhile, the sampling efficiency of point set diffusion is superior. The overall presentations of both the methodology and experiments are excellent.
weaknesses: In my opinion, the main weakness, or the most improvement-needed part of the paper, lies in the modeling and experiments of ordered point processes:

1. An important characteristic of the ordered point processes (TPPs or STPPs) is the dependence between future events and past events, which is not considered in the model. The proposed method seems to only consider the first-order statistics of the data (the event intensity/density), and treat these statistics at certain times or locations as fixed values to be learned by the model. For example, if the training data set contains multiple event trajectories sampled over the horizon of $[0, T]$, then the model will assume $p(T/2)$ (the event density at $T/2$) is fixed and to be learned. However, in TPPs, the $p(T/2)$ depends on the history (observation before $T/2$), and is different in each realization of the event trajectory, which violates the assumption of the diffusion model.

2. Although a conditional sampling method is proposed in the paper (Algorithm 3.1), I am wondering about its effectiveness in practice. First, what the $q(X_{t-1}|X_{0}^{c})$ is (line 287) remains unknown. Meanwhile, the (technical/practical) reason for using this conditioning is not shown. The results in Figure 7 are not convincing enough. To me, even though the authors claim that they are solving conditioning tasks and are visualizing the predicted densities for events from different trajectories (panels at the bottom), these density plots would look similar if we overlap them with each other. In other words, I think the model only predicts an averaged event intensity over space, and it has little connection with the conditioned samples.

3. An alternative to prove the effectiveness of conditional generation is to show the predicted intensity/density of events at different times, given a trajectory from the pinwheel dataset. This is the same idea as Figure 5 in [1]. The difference between density functions at different times is more significant and would help validate the conditional sampling method.

4. The conditional sampling task is only experimented with in the spatial domain. An example of showing the evolution of the predicted conditional density of a pinwheel trajectory can support the claim of effective conditional sampling in an ordered (temporal) domain.

5. I am also concerned that there is no log-likelihood metric reported in the paper. The metrics used in the paper are about the first-order characteristics of the data, on which I believe the proposed Point Set Diffusion can perform well. However, they cannot fully reflect the model's fit to the data when second-order data dependencies are involved (e.g., in ordered point processes). On the other hand, the log-likliehood is still the golden standard to suggest the model's goodness-of-fit to the data when it comes to conditional models or tasks [2][3]. Other point process studies that use the diffusion model will also report the data log-likelihood when evaluating the model [4][5]. I am curious about the proposed model's performance on the log-likelihood metric.

Again, I acknowledge and respect the authors' contribution to the proposed method, and I hope the above questions can be properly answered or addressed.

---
[1] Chen, Ricky TQ, Brandon Amos, and Maximilian Nickel. "Neural Spatio-Temporal Point Processes." International Conference on Learning Representations.

[2] Daryl J Daley, David Vere-Jones, et al. An introduction to the theory of point processes: volume I: elementary theory and methods. Springer, 2003

[3] Reinhart, Alex. "A review of self-exciting spatio-temporal point processes and their applications." Statistical Science 33.3 (2018): 299-318.

[4] Dong, Zheng, Zekai Fan, and Shixiang Zhu. "Conditional Generative Modeling for High-dimensional Marked Temporal Point Processes." arXiv preprint arXiv:2305.12569 (2023).

[5] Yuan, Yuan, et al. "Spatio-temporal diffusion point processes." Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining. 2023.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: 1. The sampling time and quality of the diffusion model are directly related to the number of forward/backward steps, which I could not find in the paper. Could the authors provide some ablation study on the number of steps, e.g., how the sampling time and quality grow with it?
soundness: 3
strengths: 1. The paper is overall well-written and easy to follow. The basic concepts are introduced clearly with consistent notations. The forward process, backward process, and final sampling algorithms are well explained. Illustrations (Figure 1-3) are very clear for readers to follow the workflow of the proposed Point Set Diffusion model. The datasets and metrics are also clear in the experiment section. 

2. The idea of leveraging diffusion models to generate the whole point process is intriguing, and it is quite different from the common approaches that use autoregressive models with parameterized intensity functions that suffer from sampling speed and are restricted to forecasting tasks. Numerical results are very promising to support the efficacy of the proposed model.
summary: The paper proposes a Point Set Diffusion model for conditioned and unconditioned generations of point processes (spatial, temporal, and spacial-temporal) without intensity functions. The model treats the latent space of the point process as a whole and applies diffusion to learn how to generate point processes from noise (unconditioned) and conditioning masks (conditioned). At the training phase, the point process is passed through a forward process that gradually thins the original points and adds points from a noise point process. Then, a parameterized model is trained for the backward process that gradually predicts the points in the last timestep conditioned on the current timestep and thins the noise points in the current point process. After training, both conditioned and unconditioned sampling procedures are provided. Numerical experiments illustrate that the proposed Point Set Diffusion model achieves much faster sampling speed than intensity-based autoregressive models. Moreover, it outperforms several baseline autoregressive models on various SPP and STPP tasks, especially density estimation tasks.
weaknesses: 1. Currently there are very few baseline algorithms, e.g., for SPP conditional generation there is only one baseline, and for STPP forecasting there are only two. It would be more convincing to compare with more baseline models, or to provide more evidence that the current baselines are already SOTA (which I believe they are).

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. In the experiments, how are $\alpha_t$, $\beta_t$, and $T$ set?

2. The proposed model generalizes the Add-Thin model to general metric spaces. When modeling univariate TPPs, how does the performance of the proposed model compare to Add-Thin in both unconditional and conditional sampling scenarios?
soundness: 3
strengths: 1. This paper generalizes the Add-Thin model to define a model for point processes on general metric spaces, enhancing the model's applicability and promising future prospects.
2. The idea is sound and well-founded. The paper is overall well-written and easy to follow.
3. Experiments show that the proposed model achieves state-of-the-art results on both conditional and unconditional tasks while enabling faster sampling.
summary: This paper proposes a diffusion-based latent variable model for general point processes on metric spaces. By learning stochastic interpolations between data and noise point sets, the model enables efficient, parallel sampling and flexible generation for complex tasks on the metric space. Experiments on synthetic and real-world datasets show that the proposed model achieves state-of-the-art results in unconditional and conditional tasks for spatial and spatio-temporal point processes.
weaknesses: 1. It would be helpful to discuss the connections between the proposed model and the Add-Thin model when modeling univariate temporal point processes.
2. In the conditional sampling, the definition of $q(X_{t-1} | X_{0}^c)$ in line 287 was not provided.

Typo: $X_{t+1}^{\text{thin}}$ and $X_{t}^{\text{thin}}$ in Eq.(9) should be $X_{t+1}^{\varepsilon}$ and $X_{t}^{\varepsilon}$.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - Can you add a discussion on learning time? How does the computational complexity increase, especially with more data points?
- How can hyperparameters (e.g. number of diffusion steps T or noise scheduling) be determined?
- Please tell me more about the limitation of the proposed method. For example, how robust is the proposed method in situations where there is little data? Also, will the interpretability of the proposed method be lower than parametric methods (e.g., DNN-based Hawkes processes), or will the number of sensitive hyperparameters increase by using diffusion models as a base?
soundness: 3
strengths: - The approach of modeling point processes using the diffusion model is interesting.
- Efficient sampling is achieved by making effective use of thinning.
- The effectiveness of the proposed method is evaluated on artificial and real data.
- The manuscript is well-written.
summary: This paper proposes a novel diffusion-based approach to model point processes without relying on traditional intensity functions. This model is characterized by its ability to efficiently and flexibly generate point sets through stochastic interpolation between data and noise sets. Experiments on synthetic and real-world datasets demonstrate that the model achieves state-of-the-art performance in generating spatial and spatiotemporal point processes, significantly outperforming existing methods in terms of speed and accuracy of sample generation.
weaknesses: - There is not enough discussion about computational complexity.
- Not very clear on how to set hyperparameters.
- No mention of the effectiveness of the method with respect to the amount of data.
- No discussion of limitation.

Meta-review: This paper proposes an interesting diffusion-based method for modeling and sampling point sets of point process distributions. For spatio-temporal point processes (STPP) the dominant approach is to apply autoregressive models that move points across time. This work compares to those methods, showing that they can sample point sets much more efficiently by parallelizing across time. Expert reviewers felt that the work was interesting, taking a new approach than prior work, and showing some clear advantages experimentally. We would be happy for this paper to be presented at ICLR. The authors already updated their appendix to address questions raised by the reviewers. We would encourage them to continue refining the work before the final version of the paper, and adding further experiments that could not be run due to time limitations.

Decision: Accept (Poster)