Particle gradient descent model for point process generation

Abstract

This paper presents a statistical model for stationary ergodic point processes, estimated from a single realization observed in a square window. With existing approaches in stochastic geometry, it is very difficult to model processes with complex geometries formed by a large number of particles. Inspired by recent works on gradient descent algorithms for sampling maximum-entropy models, we describe a model that allows for fast sampling of new configurations reproducing the statistics of the given observation. Starting from an initial random configuration, its particles are moved according to the gradient of an energy, in order to match a set of prescribed moments (functionals). Our moments are defined via a phase harmonic operator on the wavelet transform of point patterns. They allow one to capture multi-scale interactions between the particles, while controlling explicitly the number of moments by the scales of the structures to model. We present numerical experiments on point processes with various geometric structures, and assess the quality of the model by spectral and topological data analysis.

Appendices

Appendix A: Proof of Theorem 1

In order to prove Theorem 1, we need to formally define Eq. (9). Recall that in this section and in what follows, $W_s$ is interpreted as endowed with the addition and scalar multiplication modulo $W_s$.

For $\mu \in \mathbb M^s$ and any $x \in \text {Supp}(\mu )$, we define the following functions:

$$\begin{aligned}&\begin{array}{ccccc} {h_{x}^{\mu }}&{} : &{} \mathbb R^2 &{} \longrightarrow &{} \mathbb M^s \\ &{} &{} y &{} \longmapsto &{} \displaystyle {\mu - \delta _x + \delta _{x+y},} \\ \end{array}\\&\begin{array}{ccccc} K_{x}^{\mu }&{} : &{} \mathbb R^2 &{} \longrightarrow &{} \mathbb C^d \equiv \mathbb R^{2d} \\ &{} &{} y &{} \longmapsto &{} \displaystyle {K\circ h_{x}^{\mu }(y)}, \\ \end{array}\\&\begin{array}{ccccc} E_{x}^{\mu }&{} : &{} \mathbb R^2 &{} \longrightarrow &{} \mathbb R^+ \\ &{} &{} y &{} \longmapsto &{} \displaystyle {E_{\bar{\phi }}\circ h_{x}^{\mu }(y)}. \end{array} \end{aligned}$$

The function $K_{x}^{\mu }$ can be complex valued. However, as our energy function is the square Euclidean norm, it is equivalent to consider that $K_{x}^{\mu }$ has values in $\mathbb R^{2d}$. Moreover, we assume in what follows that the function K is such that for all $\mu \in \mathbb M^s$ and all $x \in \text {Supp}(\mu )$, $K_x^\mu$ is differentiable. We can then define from chain rule, for any $\mu \in \mathbb M^s$ and any $x \in W_s$

$$\begin{aligned} \nabla _xK(\mu ) := \left\{ \begin{array}{ll} Jac[K_x^\mu ](0) &{} \text{ if } x \in \text {Supp}(\mu )\\ 0 &{} \text{ otherwise } \end{array} \right. , \end{aligned}$$

(18)

and

$$\begin{aligned} \nabla _xE_{\bar{\phi }}(\mu ) := \left\{ \begin{array}{ll} Jac[E_x^\mu ](0) &{} \text{ if } x \in \text {Supp}(\mu )\\ 0 &{} \text{ otherwise } \end{array} \right. , \end{aligned}$$

(19)

where Jac[f] denotes the Jacobian matrix of the function f. When $x \in \text {Supp}(\mu )$, the chain-rule gives $Jac[E_x^\mu ](0) = (\nabla _xK(\mu ))^t(K(\mu )-K(\bar{\phi }))$. We can now give the proof of Theorem 1.

Proof

We are going to show that if $\Phi _n$ follows a distribution invariant to T, then $\Phi _{n+1} = G_{\bar{\phi }} (\Phi _n)$ also follows a distribution that is invariant to T. The gradient descent procedure thus produces a sequences of measures $\Phi _n$ that are all invariant to T because the initial random measure $\Phi _0$ is invariant to T.

Denote $G_{\bar{\phi }} (\mu )$ by the measure configuration transported from $\mu$, by performing one gradient-descent step on the energy $E_{\bar{\phi }}$. More precisely, for $\mu = \sum _i \delta _{x_i}$, we define for a fixed $\gamma >0$, the gradient-descent step by

$$\begin{aligned} G_{\bar{\phi }} (\mu ) := \sum _i \delta _{x_i - \gamma \nabla _{x_i} E_{\bar{\phi }} (\mu ) } . \end{aligned}$$

For any transform $T x = Ax + b$ on $W_s$, where A is an orthogonal matrix A with entries in $\{-1,0,1\}$, and $b\in W_s$. As A is a linear transformation on the torus $W_s$, $\forall x , y \in W_s$, $A(x+y) = Ax + Ay$. We shall first prove that

$$\begin{aligned} G_{\bar{\phi }} (T_{\#} \mu ) = T_{\#} G_{\bar{\phi }} ( \mu ). \end{aligned}$$

(20)

Let $y_i = T x_i$, then by definition,

$$\begin{aligned} G_{\bar{\phi }} ( T_{\#}\mu ) = \sum _i \delta _{ y_i - \gamma \nabla _{y_i} E_{{\bar{\mu }}} (T_{\#} \mu ) } , \end{aligned}$$

and

$$\begin{aligned} T_{\#} G_{\bar{\phi }} ( \mu ) = \sum _i \delta _{T( x_i - \gamma \nabla _{x_i} E_{\bar{\phi }} (\mu ) ) }. \end{aligned}$$

We are going to show that for each i-th particle, $y_i - \gamma \nabla _{y_i} E_{\bar{\phi }} (T_{\#} \mu ) = T( x_i - \gamma \nabla _{x_i} E_{\bar{\phi }} (\mu ) )$. This implies that (20) is correct. The key is to show that

$$\begin{aligned} A \nabla _{x_i} E_{\bar{\phi }} (\mu ) = \nabla _{y_i} E_{\bar{\phi }} (T_{\#} \mu ) \end{aligned}$$

(21)

which will imply that $\forall i$,

$$\begin{aligned}&T( x_i - \gamma \nabla _{x_i} E_{\bar{\phi }} (\mu ) ) = A x_i - \gamma A \nabla _{x_i} ,\\&E_{\bar{\phi }} (\mu ) + b = y_i - \gamma A \nabla _{x_i} E_{\bar{\phi }} (\mu ) = y_i - \gamma \nabla _{y_i} E_{\bar{\phi }} (T_{\#} \mu ) . \end{aligned}$$

To show (21), we recall that by the definitions in Sect. 3.1,

$$\begin{aligned}&\nabla _{x_i} E_{\bar{\phi }} (\mu ) = Jac(K_{x_i}^\mu )(0)^t ( K ( \mu ) - K(\bar{\phi }) ) , \end{aligned}$$

(22)

$$\begin{aligned}&\nabla _{y_i} E_{\bar{\phi }} (T_{\#} \mu ) = Jac(K_{y_i}^{T_{\#} \mu })(0)^t ( K (T_{\#} \mu ) - K(\bar{\phi }) ) , \end{aligned}$$

(23)

with $Jac(K_{y_i}^{T_{\#} \mu })(0) = Jac(K_{y_i}^{T_{\#} \mu } \circ A \circ A^{-1})(0) = Jac(K\circ h_{y_i}^{T_{\#}\mu } \circ A)(0)A^{-1}.$ Furthermore, $\forall x \in W_s$,

$$\begin{aligned}&K\circ h_{y_i}^{T_{\#}\mu } \circ A (x) \nonumber \\ {}&= K(T_{\#}\mu - \delta _{y_i} + \delta _{y_i + Ax}) = K(T_{\#}(\mu - \delta _{x_i} + \delta _{x_i + x})) \nonumber \\&= K(\mu - \delta _{x_i} + \delta _{x_i + x}) = K \circ h_{x_i}^\mu (x), \end{aligned}$$

(24)

where we used the fact that T is affine, and the invariance of K w.r.t. T. The equality in (21) follows directly from (22), (23), (24) and the fact that $A^{-1} = A^t$. From (21), we conclude that (20) holds.

Based on (20), it remains to show that for any Borel set $\Gamma$ on $W_s$, $P ( \Phi _n \in T_{\#}^{-1} \Gamma ) = P ( \Phi _n \in \Gamma )$. This can be shown by induction, since by assumption it holds at $n=0$: assume now that this statement holds at $n \ge 0$, then we have,

$$\begin{aligned} P ( \Phi _{n+1} \in T_{\#}^{-1} \Gamma )&= P ( G_{\bar{\phi }} (\Phi _n) \in T_{\#}^{-1} \Gamma )\nonumber \\&= P ( T_{\# } G_{\bar{\phi }} (\Phi _n) \in \Gamma ) \nonumber \\&= P( G_{\bar{\phi }} ( T_{\# } \Phi _n) \in \Gamma ) \nonumber \\&= P( G_{\bar{\phi }} ( \Phi _n) \in \Gamma )\end{aligned}$$

(25)

$$\begin{aligned}&= P ( \Phi _{n+1} \in \Gamma ) . \end{aligned}$$

(26)

The second last equality in (26) is due to the invariance of $\Phi _n$, i.e.

$$\begin{aligned} P( G_{\bar{\phi }} ( T_{\# } \Phi _n) \in \Gamma )&= P ( \Phi _n \in T_{\#}^{-1} G_{\bar{\phi }}^{-1} \Gamma )\\&= P ( \Phi _n \in G_{\bar{\phi }}^{-1} \Gamma )\\&= P( G_{\bar{\phi }} ( \Phi _n) \in \Gamma ). \end{aligned}$$

$\square$

Appendix B: Fourier spectrum and power spectrum

We define the discrete Fourier transform (DFT) $F_m(\mu )$ of a counting measure $\mu =\sum _u\delta _{x_u}\in \mathbb M^s$ on the (square) window $[-s,s[^2$ at integer frequency $m \in \mathbb Z^2$ by

$$\begin{aligned} F_m(\mu ):=\int _{ W_s} e^{-i\pi mx/s}\,\mu (dx)=\sum _{u} e^{-i\pi m x_u/s}. \end{aligned}$$

Observe, $F_m(\mu )$ at frequency $m=(0,0)$ specifies the number of points of the measure $\mu$ on $W_s$. The empirical Fourier spectrum (or power spectrum ) is often defined by taking the square modulus of the Fourier coefficients $F_m(\mu )$; $U_m(\mu ):=|F_m(\mu )|^2$. Note that $|F_m(\mu )|^2$, and consequently $U_m(\mu )$ is invariant with respect to (circular) translations of $\mu$ on $W_s$. By selecting the frequencies in a limited range $m \in \Gamma _F \subset \mathbb Z^2$, one obtains a translation-invariant Fourier spectrum. As we shall focus on isotropic point processes, we further reduce the variance of our statistics by averaging Fourier coefficients along frequency orientations. More precisely, let us define ${\tilde{\Gamma }}_F := \{ \lfloor |m| \rfloor , \, m \in \Gamma _F \}$. For each $k\in \tilde{\Gamma }_F$, we define ${\tilde{U}}_k(\mu ) := \frac{1}{\#k}\sum _{\begin{array}{c} m \in \Gamma _F \\ \lfloor |m| \rfloor = k \end{array}} U_m(\mu )$, where $\#k$ denotes the cardinal of $\{ m \in \Gamma _F : \lfloor |m| \rfloor = k \}$. The radial power spectrum P(k) is the expectation of ${\tilde{U}}_k(\mu )$ for $k \in N = \{ 1,2,3,\ldots \}$ when $\mu$ follows some distribution, divided by the intensity of the process (estimated over 10 realizations).

Appendix C: Relaxing the assumptions on the data

In this paper, in order to present our model in a simple setting, strong theoretical assumptions have been made on the data. However, in real world applications, the data will most likely not satisfy these assumptions. This sections presents ideas on how to adapt our model in such cases.

Non-periodic boundaries Recall that our descriptor, defined in (14), applies periodic boundary correction to point patterns in a square window. If the structure of the observed pattern is not periodic, one can modify the descriptor by applying non-periodic integrals in (14) over some smaller window. In particular, we suggest a scale-dependent reduction of the integration window, pertinent when the wavelet $\psi$ has a compact (or approximately compact) spatial support. Specifically, we consider a new descriptor ${\tilde{K}}$ by considering the integrals in (14) with $i=(\lambda ,k,\lambda ',k', \tau ') \in \Gamma _H$ over smaller windows $W_{s_i}\subset W_s$, such that boundary effects are negligible. Our current software can also handle such non-periodic boundary conditions.

More general observation windows In this paper, we considered that the observed pattern lies in a square observation window. If this is not the case, one could use a similar idea to the non-periodic case: embed the observation window in a square window and considering integrals in (14) over the observation window.

Non stationary process In Koňasová and Dvořák (2021), the authors focus on building a model for non stationary point processes inspired by Tscheschel and Stoyan (2006). Similarly, one might adapt our method to model non stationary processes. This could be done by modifying two aspects of the method. First, the initial distribution $\Phi _0$ (cf. Sect. 3.1) could be chosen as a non stationary Poisson point process, estimating the intensity with a kernel estimator, such as in Koňasová and Dvořák (2021). In addition, as pointed out in Koňasová and Dvořák (2021), the descriptor should be adapted not to be translation invariant. This could be done, for example, by applying local integrals over patches of the observation window in (14) (this requires some notion of “local stationarity” of the process). Another method that may be useful in this scenario is the regularization proposed in Brochard et al. (2020), where a regularization term is added to the energy. This term consists of the (Sliced Wasserstein) distance (Rabin et al. 2011) between the initial configuration and the current configuration (the one being optimized). By adding this regularization term to the energy, the points of the configuration are forced not to move too far away from the initial configuration, which could help preserve the non stationarity of the initial distribution in the distribution of the model.

Table 6 Discussion of the important parameters of our model

Processes in other dimensions While we focus in this paper on planar point processes, our approach can readily be extended to any dimensions. To model point processes in other dimensions such as 1d or 3d, one can consider similar type of wavelets proposed in the literature (Chenouard and Unser 2011; Brumwell et al. 2018).

Appendix D: List of important parameters of our model

In Table 6, we discuss the main parameters of our model in three categories. The first two categories are the parameters that are relatively standard to consider in most existing methods such as Tscheschel and Stoyan (2006). The third category is more specific to our model, which involves the discretization step, and the final blurring step.