Is Score Matching Suitable for Estimating Point Processes?

Thank you so much for reviewing our paper. We answer your questions below.

The argument to choose their weighting function $h$ is definitely better than simply picking one. That being said, it still not clear how close to optimal their choice of h is given they are only optimizing one term in the objective $C_{\bf h}$. The paper would be a bit stronger if they could include a different (natural) weighting function to show their new choice is an improvement.

We thank the reviewer for pointing this out. To address your concern, we carry out experiments and compare three different weight functions that satisfy the criterion (Eq. 15) for AWSM on Hawkes process. Other than the near-optimal weight function $\bf{h}^0$ introduced in our paper, we also consider a natural weight function $\bf{h}^1$ defined as,

$h_n^1(t_n)=(t_n-t_{n-1})(T-t_{n}).$

We also consider another valid weight function: the square root of $\bf h^1$ denoted as $\bf h^2$,

$h_n^2(t_n)=\sqrt{(t_n-t_{n-1})(T-t_n)}.$

All three weight functions can be applied in AWSM to recover ground-truth parameters, however with different convergence rates. We carry out experiments on synthetic data for exponential-decay model with the same setting as section 6.2 in our paper. We measure their MAE for different sample sizes in Figure 1 of the attached PDF and find that $\bf h^0$ does achieve the best results among the three weight functions.

We hope these additional experimental results can enhance the credibility of our paper and address your concerns. We would definitely add these results in the camera ready.

Given the choice of $\bf h$ optimizes $C_{\bf h}$ but ignores $\gamma(\bf h,A,B)$, how sensitive is $\Gamma$ to the choice of $\bf h$? This could be demonstrated empirically or theoretically, but it would be good to have some analysis of this to argue why this weighting function is decent.

We thank the reviewer for pointing this out. This is indeed a tough question that we do not yet have a satisfying answer. For specific parametric models, $\dot A_n(\mathcal T), \dot B_n(\mathcal T)$ can be computed analytically (see line 478-480) and then $\Gamma(\mathbf{h},A,B)$ can be computed via Monte Carlo. Then we can study how sensitive $\Gamma(\mathbf{h},A,B)$ is to $\bf h$. For general models, especially when $\psi_{\theta}$ is a deep neural network like THP or SAHP, $\Gamma(\mathbf{h},A,B)$ is intractable to compute.

However, heuristically speaking, our choice of near-optimal weight function $\bf h^0$ should be a good choice even concerning $\Gamma(\mathbf{h},A,B)$. To make $\Gamma(\mathbf{h},A,B)$ small, a natural idea is to make $|h_n(\mathcal T)|$ and $|\frac{\partial}{\partial t_n}h_n(\mathcal T)|$ small. The weight function we chose and its derivative have relatively positive low powers with respect to $t_n$, therefore making $|h_n^0(\mathcal T)|$ and $|\frac{\partial}{\partial t_n}h_n^0(\mathcal T)|$ small. For weight functions $h_n^1(\mathcal T)=(T-t_n)(t_n-t_{n-1})$, the power w.r.t. $t_n$ is two. And for $h^2_n(\mathcal T)=\sqrt{(T-t_n)(t_n-t_{n-1})}$, its derivative is $\frac{\partial}{\partial t_n}h_n^2(\mathcal T)=\frac{1}{2}\frac{T-t_n-(t_n-t_{n-1})}{\sqrt{(T-t_n)(t_n-t_{n-1})}}$ , the numerator is usually a bounded quantity and the denominator may be close to zero, making it large. In conclusion, $\bf h^0$ is a better choice compared with $\bf h^1$ or $\bf h^2$ concerning $\Gamma(\mathbf{h},A,B)$.

Equation 1: Can we at least have a sentence explaining the paper will be slightly abusing notation by referring to $p(\mathcal T)$ as the conditional distribution of p given $N_T=|T|$?

We thank the reviewer for pointing this out. Here since $N_T$ is random, by the notation $p(\mathcal T)$, we are referring to the probability density or the likelihood of observing $N_T$ events $t_1,\ldots, t_{N_T}$ in $[0,T]$. For Poisson process, the conditional distribution of event time $t_1,\ldots, t_{N}$ given $N_T = N$ is,

$p(t_1,\ldots, t_N|N_T = N) = \frac{\prod_{n=1}^N \lambda(t_n)}{\int_{0\leq t_1\dots\leq t_N\leq T}\prod_{n=1}^{N}\lambda(t_n)dt_1\ldots dt_N}, 0\leq t_1\ldots\leq t_{N}\leq T.$

We acknowledge that our paper may not employ the rigorous notation typically used in the fields of probability and stochastic processes. Instead, we have opted for notations that we believe are accessible to readers from both the computer science and statistics communities. We deeply apologize for any confusion this may have caused.

Thank you so much for reviewing our paper. We answer your questions below.

Denoising score matching for point processes is not considered.

We thank the reviewer for pointing this out. Indeed, denoising score matching (DSM) is not considered in our paper because we mainly focus on correcting the original score matching. However, to address your concerns, we have carried out experiments on DSM. We can conclude that DSM is inferior to AWSM in terms of both efficiency and accuracy.

We deploy DSM on THP and SAHP. We use the DSM loss as in [1]. For observed timestamps $t_n^{(m)}$ in $m$-th sequence, we sample $L$ noise samples $\tilde t_{n,l}^{(m)}=t_n^{(m)}+ \epsilon_{n,l}^{(m)}, l = 1,\ldots, L,$ where $\text{Var}(\varepsilon_{n,L}^{(m)})=\sigma^2$ and get the DSM objective:

$\hat {\mathcal J}(\theta)=\frac{1}{M}\sum_{m=1}^M\sum_{n=1}^{N_m}\sum_{l=1}^{L}\frac{1}{2L}[\psi_\theta(\tilde t_{n,l}^{(m)})+\frac{\varepsilon_{n,l}^{(m)}}{\sigma^2}]+\alpha\hat {\mathcal J}_\text{CE}(\theta),$

where $\hat {\mathcal J}_{\text{CE}}(\theta)$ is the cross-entropy loss as line 191-192 in our paper. We carried out experiments on synthetic data where $T$ is known. We compared the TLL, ACC, and running time of DSM with MLE and AWSM. The comparison results are shown in Table 1 in the attached PDF. DSM performs the worst among the three estimators, and AWSM is the fastest.

The reasons why DSM performs poorly are as follows:

1. In terms of accuracy, as discussed in Section 5 of [2], DSM is biased when $\sigma > 0$, while our AWSM is unbiased and produces consistent estimates. Therefore, it is not surprising to see that AWSM achieves better TLL and ACC than DSM.

2. In terms of efficiency, previous work states that DSM is faster than SM because the Hessian matrix is expensive to compute, while DSM avoids this. However, as presented in [3], ASM (or AWSM) also avoids computing the Hessian matrix by replacing it with a one-dimensional partial derivative $\frac{\partial}{\partial t_n}\psi_{\theta}(t_n|\mathcal F_{t_{n-1}})$. Therefore, DSM is not more efficient than AWSM for avoiding the computation of the Hessian matrix. On the contrary, DSM can be slower than AWSM because DSM requires $L$ noise samples for each timestamp and must compute their scores.

For Hawkes processes, the problem can be solved by applying log-normalization to transform the bounded temporal domain into an unbounded one [1].

This is a misunderstanding. Apply log-normalization does not provide any solution to the problem, we can prove that the integration by parts still produces an intractable term after taking this transformation. Therefore, using ASM on a log-transformed sequence still results in wrong estimation. If one hopes to get consistent estimation after log transformation, a suitable weight function is still needed. This is summarised as follows.

We illustrate this using the simplest example both theoretically and empirically. Consider a homogeneous Poisson process with intensity $\lambda ^*$ with observations $t_1, \ldots, t_{N_T}$ over $0\leq t_1\leq \ldots \leq t_{N_T}\leq T$. Now we consider

1. First apply log transformation to timestamps then using ASM. The estimate is denoted as $\hat \lambda_{\text{ASM,Log}}$.

2. First apply log transformation to intervals then using ASM on intervals. The estimate is denoted as $\hat \lambda_{\text{ASM,Log Interval}}$.

For log-transformation on timestamp, we have $x_n=\log t_n$. Therefore the conditional pdf of $x_n$ given $\mathcal F_{x_{n-1}}$ is:

$p(x_n|\mathcal F_{x_{n-1}})=\lambda \exp\left[-\lambda \exp(x_n)\right]\exp(x_n),$

and we can compute conditional score accordingly. We plug these terms into the implicit ASM objective and get the estimate. Suppose we observed $M$ sequences of $t_1^{(m)}, \ldots, t_{N_m}^{(m)}$, after derivation we get an analytical form of $\hat \lambda_{\text{ASM,Log}}$:

$\hat \lambda_{\text{ASM,Log}} =2\frac{\sum_{m=1}^M\sum_{n=1}^{N_m}\exp(x_n)}{\sum_{m=1}^M\sum_{n=1}^{N_m}\exp(2x_n)}.$

Similarly, we consider transformation on intervals $\tau_n:= t_n-t_{n-1}, n\geq 2, \tau_1:=t_1, y_n = \log (\tau_n)$ and consider the score $\psi(y_n|\mathcal F_{y_{n-1}})$ then apply ASM. The estimate in this case is:

$\hat \lambda_{\text{ASM,Log Interval}} =2\frac{\sum_{m=1}^M\sum_{n=1}^{N_m}\exp (y_n)}{\sum_{m=1}^M\sum_{n=1}^{N_m}\exp(2y_n)}.$

Both estimates are wrong and can never recover the true parameter regardless of sample size. To compare, after transformation, suitable weight functions can be added. For log transformation on timestamps, we use weight function $h_n(x_n)=(x_n-x_{n-1})(\log T-x_n)$. For log transformation on intervals, we use $h_n(y_n) = \log [T-t_{n-1}]-y_n$. The corresponding estimates are $\hat \lambda_{\text{AWSM,Log}}$ and $\hat \lambda_{\text{AWSM,Log Interval}}$.

We measure the MSE of these four estimates. Results are shown in Table 2 in the attached PDF. It is easy to see that as sample size increases, MSE of $\hat \lambda_{\text{ASM,Log}}$ and $\hat \lambda_{\text{ASM,Log Interval}}$ remains large and unchanged, showing that these two estimators are wrong regardless of sample size.

Including more metrics would be helpful.

We will consider adding MSE to our revised paper.

Did the author observe any instability during training?

For our weighted score matching, we did not notice instability.

[1] SMURF-THP: score matching-based uncertainty quantification for transformer Hawkes process

[2] A connection between score matching and denoising autoencoders

[3] Autoregressive score matching

We greatly thank the reviewer for your constructive feedback.

For the log-normalization, could you kindly provide the detailed derivation for obtaining $\hat \lambda_{\text{ASM,Log}}$ and $\hat \lambda _{\text{ASM,Log Interval}}$?

Sure, we are more than delighted to discuss this with you.

For log-normalization on timestamps, since the transformed variable $x_n$ has conditional pdf,

$p(x_n|\mathcal F_{x_{n-1}})=\lambda \exp\left[-\lambda \exp(x_n)\right]\exp(x_n).$

The conditional score and its partial derivative w.r.t. $x_n$ will be,

$\psi(x_n|\mathcal F_{x_{n-1}}) = -\lambda \exp(x_n) + 1,$

$\frac{\partial}{\partial x_n}\psi_{\theta}(x_n|\mathcal F_{x_{n-1}})=-\lambda \exp(x_n).$

Now we perform ASM on the transformed variable $x_n$. Denote $\mathcal X = (x_1,\ldots, x_{N_T})^T$ as the vector of transformed timestamps, the implicit ASM objective (Eq. 13) for $\mathcal X$ will be,

$\mathcal J_\text{ASM}(\theta)=\mathbb E_{p(\mathcal X)}[\sum_{n=1}^{N_T}\frac{1}{2}\psi^2_\theta(x_n|\mathcal F_{x_{n-1}})+\frac{\partial}{\partial x_n}\psi_\theta(x_n|\mathcal F_{x_{n-1}})].$

We plug in the conditional score we just computed into the above equation and get,

$\mathcal J_{\text{ASM}}(\theta) =\mathbb E_{p(\mathcal X)}[\sum_{n=1}^{N_T}\frac{1}{2}\exp(2x_n)\lambda^2-\sum_{n=1}^{N_T}2\exp(x_n)\lambda ]+C,$

where $C$ is a constant not containing $\theta$.

The empirical objective based on samples $t_1^{(m)},\ldots t_{N_m}^{(m)}, m=1,\ldots M$ will be,

$\hat J_\text{ASM}(\theta)=\frac{1}{M}\sum_{m=1}^M\sum_{n=1}^{N_m}[\frac{1}{2}\exp(2x_n^{(m)})\lambda^2-2\exp(x_n^{(m)})\lambda].$

This is quadratic with respect to $\lambda$, and its minimal $\hat \lambda_{\text{ASM,Log}}$ can be computed analytically as $\hat \lambda_{\text{ASM,Log}} = 2\frac{\sum_{m=1}^M\sum_{n=1}^{N_m}\exp(x_n^{(m)})}{\sum_{m=1}^M\sum_{n=1}^{N_m}\exp(2x_n^{(m)})}.$

For log transformation on time interval $y_n$, since the conditional pdf of $\tau_n = t_n-t_{n-1}$ is,

$p(\tau_n|\mathcal F_{\tau_{n-1}})=\lambda \exp(-\lambda \tau_n),$

thus the conditional pdf of $y_n = \log \tau_n$ will be,

$p(y_n|\mathcal F_{y_{n-1}})=\lambda \exp\left[-\lambda \exp(y_n)\right]\exp(y_n),$

it has the same analytical expression as that of $x_n$, except that $y_n$ and $x_n$ have different support. Therefore, following exactly the same steps as those of $x_n$, we can get the analytical expression of $\hat \lambda_{\text{ASM,Log, Interval}}$ in our rebuttal (with minor modification for adding superscript $(m)$).

Additionally, is it possible to express $\hat \lambda_{\text{AWSM,Log}}$ and $\hat \lambda_{\text{AWSM, Log Interval}}$ in closed form?

Yes, it is possible. We take $\hat \lambda_{\text{AWSM, Log}}$ as an example. The derivation of $\hat \lambda_{\text{AWSM, Log Interval}}$ is basically the same.

First, the implicit AWSM objective as Eq. 17 in our paper for $\mathcal X$ will be,

$\mathcal J_\text{AWSM}(\theta)=\mathbb E_{p(\mathcal X)}[\sum_{n=1}^{N_T}\frac{1}{2}\psi^2_\theta(x_n|\mathcal F_{x_{n-1}})h_n(\mathcal X)+\frac{\partial}{\partial x_n}\psi_\theta(x_n|\mathcal F_{x_{n-1}})h_n(\mathcal X)+\psi_\theta(x_n|\mathcal F_{x_{n-1}})\frac{\partial}{\partial x_n}h_n(\mathcal X)].$

Plug in $\psi_\theta(x_n|\mathcal F_{x_{n-1}})$ and $\frac{\partial}{\partial x_n}\psi_\theta(x_n|\mathcal F_{x_{n-1}})$ into the empirical version of the above equation and get,

$\hat J_\text{AWSM}(\theta)=\frac{1}{M}\sum_{m=1}^M\sum_{n=1}^{N_m}[\frac{1}{2}\exp(2x_n^{(m)})h_n(\mathcal X^{(m)})\lambda^2-2\exp(x_n^{(m)})h_n(\mathcal X^{(m)})\lambda - \exp(x_n^{(m)})\frac{\partial}{\partial x_n}h_n(\mathcal X^{(m)})\lambda].$

This is still quadratic w.r.t $\lambda$, and its minimal will be,

$\hat \lambda_{\text{AWSM, Log}}=\frac{\sum_{m=1}^M\sum_{n=1}^{N_m}[2\exp(x_n^{(m)})h_n(\mathcal X^{(m)})+\exp(x_n^{(m)})\frac{\partial}{\partial x_n}h_n(\mathcal X^{(m)})]}{\sum_{m=1}^M\sum_{n=1}^{N_m}\exp(2x_n^{(m)})h_n(\mathcal X^{(m)})}.$

Since we choose weight $h_n(\mathcal X)=(x_n-x_{n-1})(\log T-x_n)$ with derivative $\frac{\partial}{\partial x_n}h_n(\mathcal X)=\log T - 2x_n+x_{n-1}$. Plug them in the above equation, then we get a closed form of $\hat \lambda_{\text{AWSM, Log}}$.

Thanks so much for providing concrete answers to my queries and issues! Most of my concerns have been resolved and I have increased my score. I believe it would be valuable to include a discussion on denoising score matching, as well as some recent related work on score matching with point processes [1,2].

[1] Beyond Point Prediction: Score Matching-based Pseudolikelihood Estimation of Neural Marked Spatio-Temporal Point Process.

[2] Exploring generative neural temporal point process.

We thank greatly for the reviewer for further engaging in the discussion.

For the log interval case, each $\tau_i$ should follow the same exponential distribution, so $\hat{\lambda}_{\text{ASM,Log Interval}}$ should approximate the correct $\lambda$ with a large sample size.

This is not true. It is widely known that for an exponential distribution on $[0,\infty)$, it indeed has expectation $\frac{1}{\lambda}$ and second moments $\frac{2}{\lambda^2}$. Then the average of exponential distribution should be equal to $\frac{1}{\lambda}$ when sample is large. However, this is not the case when considering the time interval in a Poisson process on $[0,T]$.

Considering the numerator of $\hat \lambda_{\text{AWSM,Log Interval}}$, we cannot draw the conclusion that,

$\lim_{M\rightarrow \infty}\mathbb E_{p}[\frac{\sum_{m=1}^M\sum_{n=1}^{N_m}\tau_n^{(m)}}{\sum_{m=1}^MN_m}]=\frac{1}{\lambda}$

We cannot use CLT or other things here since $\sum_{m=1}^M N_m$ is random and each $\tau_n^{(m)}$ does not necessarily have the same marginal distribution since their support is different.

To illustrate this, we compute $\frac{\sum_{m=1}^M\sum_{n=1}^{N_m}\tau_n^{(m)}}{\sum_{m=1}^MN_m}$ for $M=10^{6}, T = 1$ under different $\lambda^{*}$ and show that the expectation of it is not equal to $\frac{1}{\lambda}$.

Also for the denominator, its expectation is not equal or asymptotically equal to $\frac{2}{\lambda^2}$ when taking average. Therefore, $\hat{\lambda}_{\text{ASM,Log Interval}}$ does not approximate the correct $\lambda$ with a large sample size.

We give a brief explanation as follows. For a stochastic process in a finite window $[0,T]$. Each $\tau_i$ has support $[0, T-\sum_{i\leq n}\tau_i]$ instead of $[0,\infty)$. Therefore, firstly, each $\tau_i$ has different support and does not follow the same distribution. Secondly, for each $\tau_i$ , it has support different from that of the exponential distribution, which is $[0, \infty)$.

In conclusion, when discussing a point process on a finite time window, the right bound's effect must be considered. And the case is different from sampling a fixed length time sequence over an infinite time window.

Thank you so much for reviewing our paper. We answer your questions below.

For example, there are several assumptions on the weight function in Equation 9 but it is not clear which properties are used in the proof of Proposition 3.3. For instance, Line 400 features an explicit expression for $h_1$ and $h_N$ while the expression in Equation 9 features none of these. Please clarify whether this is in fact the case. The same issues persist for the proofs of the results in Section 4.

We greatly thank you for pointing this out. Your understanding is correct, and we apologize for the typos. The proofs for Propositions 3.3 and 4.3 are from our previous edition, where we specified one weight function instead of all weight functions that satisfy the condition. We have proofread our work and ensured that the rest of our proofs remain correct.

For Propositions 3.3 and 4.3, the main parts remain correct. Only minor modifications are needed to utilize the assumptions stated in the main text, ensuring that the proofs hold for general weight functions. We have corrected the necessary parts as follows.

From lines 398 to 404, the corrected proof is as follows:

Using the first two equations in Eq. 9, we have

$\int [p(t_1,\ldots, t_N)\frac{\partial \log p_{\theta}(t_1,\ldots, t_N)}{\partial t_{n}}h_n(\mathcal T)]\big\vert_{t_{n}=t_{n+1}}d\mathcal T_{-n}=0, \forall n\in [N]$

$\int [p(t_1,\ldots, t_N)\frac{\partial \log p_{\theta}(t_1,\ldots, t_N)}{\partial t_n}h_n(\mathcal T)]\big\vert_{t_n=t_{n-1}}d\mathcal T_{-n}=0,\forall n\in [N].$

Therefore, the first term on the right side of the second equation of Eq. 21 will disappear, and the second term equals $-\mathbb E_{p(\mathcal T)}[\sum_{n=1}^{N_T}\frac{\partial}{\partial t_n}\psi_{\theta}(t_n)h_n(\mathcal{T})+\psi_{\theta}(t_n)\frac{\partial}{\partial t_n}h_n(\mathcal T)]$.The existence of such an expectation is due to the last two equations in Eq. 9. Therefore, we complete the proof.

We can see from the proof that, in Eq. 9, the first two equations ensure that the integration by parts trick does not produce an intractable term, and the last two equations are simply regularity conditions that ensure all terms are well-defined.

Similarly, from line 431 to 434, the corrected proof is as follows:

$\mathbb E_{p(\mathcal T)}[\sum_{n=1}^{N_T} \psi(t_n|\mathcal F_{t_{n-1}})\psi_\theta(t_n|\mathcal F_{t_{n-1}})h_n(\mathcal T)]$

$=\sum_{n=1}^\infty \int p(t_1,\ldots t_n)\psi(t_n|\mathcal F_{t_{n-1}})\psi_\theta(t_n|\mathcal F_{t_{n-1}})h_n(\mathcal T)d\mathcal T_{1:n}$

$=\sum_{n=1}^\infty \int p(\mathcal T_{:n-1})p(t_n|\mathcal F_{t_{n-1}})\psi_\theta(t_n|\mathcal F_{t_{n-1}})h_n(\mathcal T)\vert_{t_n=t_{n-1}}^{t_n=T}d\mathcal T_{:{n-1}}$

$-\sum_{n=1}^\infty \int p(t_1,\ldots, t_n)[\frac{\partial \psi_\theta(t_n|\mathcal F_{t_{n-1}})}{\partial t_n}h_n(\mathcal T) + \psi_\theta(t_n|\mathcal F_{t_{n-1}})\frac{\partial h_n(\mathcal T)}{\partial t_n}]d\mathcal T_{1:n}.$

Between the second and the third line above, we omit the steps used in the derivation of Proposition 4.1 to make it concise. For the term in the third line above, it will be eliminated using Eq. 15. For the term in the fourth line above, using Lemma B.1, we have:

$-\sum_{n=1}^\infty \int p(t_1,\ldots, t_n)\left[\frac{\partial \psi_{\theta}(t_n|\mathcal F_{t_{n-1}})}{\partial t_n}h_n(\mathcal T) + \psi_{\theta}(t_n|\mathcal F_{t_{n-1}})\frac{\partial h_n(\mathcal T)}{\partial t_n}\right]d\mathcal T_{1:n}=\\ -\mathbb E_{p(\mathcal T)}[\sum_{n=1}^{N_T}\frac{\partial \psi_{\theta}(t_n|\mathcal F_{t_{n-1}})}{\partial t_n}h_n(\mathcal T) + \psi_{\theta}(t_n|\mathcal F_{t_{n-1}})\frac{\partial h_n(\mathcal T)}{\partial t_n}].$

The existence of the expectation is ensured by the last two terms of Eq. 15.

We hope the corrected proof can address your concerns.

Thank you very much for your thoughtful review and for taking the time to carefully consider our responses. We greatly appreciate your kind words and are pleased that you find our work to be well-presented, thorough, and novel.

We will certaintly correct the proofs in the updated manuscript, as you suggested. Your feedback is valuable in improving the quality of our work, and we are grateful for your support.

Thank you once again for your constructive feedback and for increasing your rating.