K$^2$IE: Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes

We would like to thank the reviewer for the highly positive comments. Below, we provide a detailed response to each comment.

I list some possible related works for the authors' consideration, however these are definitely not essential to cite.....

We appreciate the suggestion of including these important references. "Sparse Spectral Bayesian Permanental Process with Generalized Kernel" and "Exact, Fast and Expressive Poisson Point Processes via Squared Neural Families" both propose intensity estimation methods that employ the quadratic link function, $\sigma(\cdot) = (\cdot)^2$, to ensure non-negativity of the estimators. We will appropriately cite these works in the Other Related Works section.

Is your squared loss related to equation (12) of "PSD Representations for Effective Probability Models"?

This reference proposes a quadratic kernel model of the form $f_{\text{PSD}}(x) = \sum_{n,n'=1}^N A_{nn'}k(x_n,x)k(x_{n'},x)$, which leverages the non-negativity property of the positive semi-definite matrix $A$. The model is applicable to both density estimation and intensity estimation tasks. While the reference also considers learning the model within the framework of penalized least squares loss minimization, it focuses on fitting the parameter $A$ using iterative optimization methods and does not address whether the functional form of $f_{\text{PSD}}(x)$ is optimal for that loss--indeed, it likely is not. In fact, (Marteau-Ferey et al., NeurIPS2020) showed that $f_{\text{PSD}}(x)$ minimizes a certain functional loss, $L(f(x_1),\dots, f(x_N))$, with appropriate regularization. However, this functional loss does not include the least squares loss, which involves the integral of the latent function. Although identifying the functionally optimal estimator under the least squares loss with similar regularization is an intriguing and important question, the work "PSD Representations for Effective Probability Models" does not appear essential in this specific context.

Weaknesses: As pointed out by authors in section 3.1, the method leads to intensity functions which can be negative, .... The authors use a post-hoc clipping of the intensity function below zero, after the fitting procedure, however this then clearly breaks the optimality of the solution.

In fact, applying post-hoc clipping via $\text{max}(\hat{\lambda}(x),0)$ consistently improves the accuracy of the estimator, since for any input $x$, the inequality $|\text{max}(\hat{\lambda}(x),0) - \lambda(x)| \leq |\hat{\lambda}(x) - \lambda(x)|$ holds due to the non-negativity of the true intensity function $\lambda(x)$. The reason why the clipping, despite the fact that it breaks the optimality of the solution, improves the accuracy of K$^2$IE is clear: the non-negativity condition of intensity function is not taken into consideration in the problem of minimization of the penalized least squares loss (11). For more details, see our 3rd response to Reviewer XSJm.

Question: Could you add some more discussion for equation (10), beyond providing the citations?

In response to the reviewer's suggestion, we will include the following explanation of the least squares loss (10), which we believe is satisfactory. Let $\mathbb{E}$ be the expectation regarding data points generated from the true intensity function $\lambda(x)$. Then we consider the expectation of the integrated squared loss between the estimator, $\hat{\lambda}(x)$, and the true intensity $\lambda(x)$, defined by $\mathbb{E} \Bigl[ \int_{\mathcal{X}} \bigl|\hat{\lambda}(x) - \lambda(x) \bigr|^2 dx \Bigr] = \mathbb{E} \Bigl[ \int_{\mathcal{X}} \hat{\lambda}^2(x) dx \Bigr] - 2\mathbb{E} \Bigl[ \int_{\mathcal{X}} \hat{\lambda}(x) \lambda(x) dx \Bigr] + \mathbb{E} \Bigl[ \int_{\mathcal{X}} \lambda^2(x) dx \Bigr].$ The third term can be omitted, as it is independent of the estimator. The second term can be decomposed into two parts as $2\mathbb{E} \Bigl[ \int_{\mathcal{X}} \hat{\lambda}(x) \lambda(x) dx \Bigr] = 2\mathbb{E} \Bigl[ \int_{\mathcal{X}} \hat{\lambda}(x) \sum_{n=1}^N \delta(x-x_n) dx \Bigr] + 2\mathbb{E} \Bigl[ \int_{\mathcal{X}} \hat{\lambda}(x) \bigl( \lambda(x) - \sum_{n=1}^N \delta(x-x_n) \bigr) dx \Bigr],$ where the second term vanishes due to Campbell’s theorem:

$

2\int_{\mathcal{X}} \mathbb{E} [\hat{\lambda}(x)] \lambda(x) dx - 2\sum_{n=1}^N \mathbb{E} [\hat{\lambda}(x_n)] = 2\int_{\mathcal{X}} \mathbb{E} [\hat{\lambda}(x)] \lambda(x) dx - 2\int_{\mathcal{X}} \mathbb{E} [\hat{\lambda}(x)] \lambda(x) dx = 0.

$

Putting everything together, we obtain the following (correct) interpretation of Eq. (10): $\mathbb{E} \Bigl[ \int_{\mathcal{X}} \bigl|\hat{\lambda}(x) - \lambda(x) \bigr|^2 dx \Bigr] = \mathbb{E} \Bigl[ \int_{\mathcal{X}} \hat{\lambda}^2(x) dx - 2\sum_{n=1}^N \hat{\lambda}(x_n) \Bigr] + C,$ where $C$ is the constant term.

We thank the reviewer for giving valuable comments. Below, we provide a detailed response to each question. We will include all discussions in the revised manuscript.

The equivalence between ( q(\cdot,\cdot) ) and ( h(\cdot,\cdot) ) in Equation 8 is asserted but lacks a proof ... Can this be rigorously shown without invoking path integral ...?

Yes, the equivalence between $q(\cdot,\cdot)$ and $h(\cdot,\cdot)$ can alternatively be established via Mercer's theorem, which aligns with the "rigorous" approach proposed by (Flaxman, 2017). Due to space limitations, we provide only a brief overview of the derivation here. Let $\lambda(x)$ be a function in the RKHS associated with the kernel $k(x,x')$, and let $||\lambda||^2_{H_k}$​ denote its squared RKHS norm. Then, following the notation in Eq. (7) in Section 2.2, we can express $\lambda(x) = \sum_{m=1}^{\infty} b_m \eta_m$ and $||\lambda||^2_{H_k} = \sum_m b_m^2 /\eta_m$, where $b_m$ are coefficients. Substituting this into the penalized least squares loss in Eq. (11), we obtain: $-2 \sum_{n} \lambda(x_n) + \sum_m \frac{\eta_m+1/\gamma}{\eta_m} b_m^2$. We can see that the 2nd term corresponds to the squared RKHS norm under a rescaled kernel defined by $q(x,x') = \sum_m \frac{\eta_m}{\eta_m + 1/\gamma}e_m(x) e_m(x')$. It is evident that $q(\cdot,\cdot)$ is consistent with $h(\cdot,\cdot)$ in Eq. (8). A full derivation of Theorem 1 will be included in the text.

How prevalent are negative values in practice?

Following the reviewer’s question, we conducted an analysis of how frequently K$^2$IE produces negative values using the 2D synthetic dataset $\lambda_{\text{2D}}^{1.0}$. Specifically, we evaluated the estimated intensity values at 500 x 500 grid points within the observation domain and computed the ratio of negative values. The mean ± standard deviation of this ratio across 100 trials was $0.059 \pm 0.016$. This result indicates that K$^2$IE can indeed produce negative estimates in practice—particularly in regions with sparse data—highlighting the necessity of post-hoc clipping like $\max(\hat{\lambda}(x), 0)$ in applications where negative intensity values are not permitted. It is also worth noting that when Laplace RKHS kernels are used, the equivalent kernel $h(x, x')$ is functionally non-negative in one-dimensional input settings (see Section 3.2 in (Kim, NeurIPS2024)).

What fundamentally new theoretical insight does K$^2$IE offer?

1. To the best of our knowledge, our paper is the first to prove that minimizing the least squares loss with a squared RKHS norm regularizer yields kernel intensity estimators (KIEs).

2. K$^2$IE demonstrates that the equivalent kernel used in Flaxman (2017) can serve as an edge-corrected smoothing kernel for KIE. This insight enhances computational efficiency, as it eliminates the need to solve the dual optimization problem (9) in Flaxman's model.

The claim "no model fitting" is misleading ... How does hyperparameter tuning affect computational efficiency claims?

We will rephrase "no model fitting" as "no optimization of dual coefficients". Regarding hyperparameter tuning, K$^2$IE offers a significant advantage over the reference models. Specifically, KIE and FIE require MC integration and solving a dual optimization problem for each cross-validation, respectively, whereas K$^2$IE requires neither, which is beneficial especially in multi-dimensional settings.

How does RFF rank ( M ) trade off edge-correction accuracy?

Please refer to our 2nd response to Reviewer XSJm regarding the ablation study on $M$. We want to emphasize that the edge correction (i.e., the integral over the domain in Eq. (8)) is performed exactly under the RFF approach (see the 2nd paragraph of Sec. 3.2.2). Rather, $M$ governs the approximation accuracy of shift-invariant kernels.

How does (\gamma) balance edge effects and overfitting?

We appreciate the important question. In classical KIEs, tanking $\beta \to \infty$ ($\beta^{-1}$ is the scale parameter) leads to $\hat{\lambda}(x) = \sum_n \delta(x-x_n)$. K$^2$IE exhibits the same behavior regardless of $\gamma$, but taking $\gamma \to \infty$ also yields the same solution regardless of $\beta$ in K$^2$IE. Although $\gamma$ as well as $\beta$ controls the degree of overfitting, it seems that $\gamma$ acts globally and $\beta$ locally. A thorough investigation into the distinct roles of $\gamma$ and $\beta$ remains an important direction for future work. We believe the explanation is satisfactory to the reviewer.

Experiments are limited to synthetic 1D/2D data.

Please see our 1st response to Reviewer XSJm.

Standardize notation

We will standardize the notation according to the suggestion.

How does K$^2$IE compare to modern Bayesian nonparametrics?

Please see our 3rd response to Reviewer upvZ.

does not clarify whether KIE’s edge correction was optimized

As in Sec 4, KIE's hyperparameter was optimized similarly to K2IE but based on test likelihood.

We would like to thank the reviewer for the highly positive and constructive comments, by which we are strongly encouraged. Below, we provide a detailed response to each of the comments.

Why are real-world datasets not included in the evaluation? This would be critical to support claims of practical relevance.

We focused our evaluations on synthetic datasets, which allow for precise error estimation between the true and estimated intensity functions—an appropriate setting to verify the theoretical soundness of our model. However, from the perspective of practical relevance, we fully agree with the reviewer on the importance of validating our approach on real-world datasets. In response to the reviewer’s suggestion, we have conducted an additional experiment using an open 2D real-world dataset, bei, in R package spatsta (GPL-3). It consists of locations of 3605 trees of the species Beilschmiedia pendula in a tropical rain forest (Hubbell & Foster, 1983). Following (Cronie et al., 2024), we randomly labeled the data points with independent and identically distributed marks {1, 2, 3} from a multinomial distribution with parameters $(p_1 = p_2 = 0.3, p_3 = 0.4)$, and assigned the points with label 1 and 2 to training data and test data, respectively; we repeated it ten times for evaluation. We evaluated the predictive performance of the estimators $\hat{\lambda}(x)$ based on the test least squares loss ($L_{s}$) and the negative test likelihood of counts ($L_{c}$): $L_{s}$ was computed as $L_{s} = \int_{X} \hat{\lambda}^2(x) dx - 2\sum_{n \in D_{\text{test}}} \hat{\lambda}(x_n)$, where $D_{\text{test}}$ was the test data; the observation domain $X$ was discretized into10 x 10 sub-domain {$X_1, \dots, X_{100}$}, and $L_{c}$ was computed as $L_{c} = \sum_{i=1}^{100} (\Lambda_{X_i}) - N_{X_i} \log \Lambda_{X_i} - \log(N_{X_i}!)$, where $N_{X_i}$ is the number of test data points observed in $X_i$, and $\Lambda_{X_i} = \int_{X_i} \hat{\lambda}(x) dx$. We obtained the results across 10 trials with standard errors in brackets (the lower, the better): $L_s$ = -5.74(0.39), -6.17(0.52), -5.09(0.31) for KIE, K$^2$IE, FIE; $L_s$ = 265(14.1), 278(10.8), 287(19.7) for KIE, K$^2$IE, FIE. The results show that our K$^2$IE achieved the best performance on $L_s$, but was outperformed by KIE on $L_c$, which could be because the hyperparameters were optimized based on the least squares loss and the log-likelihood for K$^2$IE and KIE, respectively. We will include the results based on 100 trials in the final manuscript.

How sensitive is the model to the number of random Fourier features (M)? An ablation study would help demonstrate robustness.

We thank the reviewer for the constructive feedback. In response to the suggestion, we conducted an ablation study on the number of random features ($2M$) using the 2D synthetic dataset $\lambda_{\text{2D}}^{1.0}$​. The integrated squared errors ($L^2$) of our K$^2$IE were 149(6.72), 74.8(12.1), 50.7(7.71), and 49.0(8.80) for $2M$ = $20$, $100$, $300$, and $500$, respectively, where standard deviations are in brackets. Similarly, the integrated absolute errors ($|L|$) were 9.86(0.25), 6.65(0.60), 5.43(0.52), and 5.31(0.57) for the same values of $M$. These results indicate that the predictive performance of K$^2$IE improves with larger $M$, and our chosen setting ($2M=500$) is sufficiently large to yield accurate estimates. Please note that, due to the limited time available during the rebuttal period, the reported results are based on 10 trials. We will include the results based on 100 trials in the final manuscript.

Can the non-negativity of the estimator be more rigorously enforced, e.g., through a post-processing projection or transformation?

As discussed in the last paragraph of Section 3.1, applying $\text{max}(\hat{\lambda}(x),0)$ or $\text{ReLU}(\hat{\lambda}(x))$ is the simplest way to enforce the non-negativity of the estimator $\hat{\lambda}(x)$. Fortunately, this operation improves the estimator's accuracy in a pointwise sense, as it satisfies $|\text{max}(\hat{\lambda}(x),0) - \lambda(x)| \leq |\hat{\lambda}(x) - \lambda(x)|$ for any input $x$ due to the fact that the true intensity function $\lambda(x)$ is always non-negative. However, when evaluating the integral $\int_{S} \lambda(x) dx$ over a compact region $S$, applying the max operation hinders closed-form integration and necessitates the use of computationally expensive Monte Carlo methods. Therefore, whether or not to employ the max⁡ operation depends on the specific application. We will incorporate the above discussion into the text. For clarity, we note that the max⁡ operation was not used in our experiments.

overlooks some recent work in Bayesian nonparametric methods...

We added a study with a scalable Bayesian model (see our 3rd response to Reviewer upvZ). We promise to cite deep kernel/GP approaches, but it would be helpful if you could point us to a few references worth citing.

We thank the reviewer for the deep understanding of our model and the constructive comments. We provide a detailed response to each comment.

Why not include a lambda>=0 constraint?, My feeling is that the paper presents a theoretically sound solution to the wrong problem.

As pointed out, one can enforce the non-negativity of the estimator by introducing lambda>=0 constraint or by applying a non-negative transformation. However, as far as we know, such constraints prevent us from obtaining efficient estimators like K$^2$IE. The main contribution of our work lies in showing that, by sacrificing strict non-negativity, one can obtain a feasible kernel-based estimator comparable to classical KIEs.

As discussed in Section 3.1, non-negativity can be enforced in a post-hoc clipping like max($\hat{\lambda}(x)$,0). However, we totally agree with the reviewer's point that the functional optimization problem (11) does not explicitly take non-negativity into account, hence not yielding an optimal "non-negative" estimator. For future work, we would like to discuss technical issues arising when non-negativity constraints are imposed.

Consider modeling the intensity as a non-negative transformation $\sigma(x)$ of a latent function $f(x)$ lying in an RKHS. Then, the functional derivative of the objective functional in Eq. (11) leads to the following equation that the optimal $\hat{f}(x)$ solves (full derivation is omitted): $\frac{1}{\gamma} f(x) + \int_{X} k(x,s)\sigma(f(s)) \sigma'(f(s)) ds = \sum_n k(x,x_n) \sigma'(f(x_n)),$ where $\sigma'(x)$ is the derivative of $\sigma(x)$. When $\sigma(x) = x$, the equation reduces to a Fredholm integral equation for which Theorem 1 provides a feasible solution. However, when $\sigma(x)$ is nonlinear, even as simple as $\sigma(x) = x^2$, deriving a feasible solution becomes non-trivial. Alternatively, one may consider enforcing non-negativity of intensity at finite virtual points {$q_1,\cdots, q_R$}, which leads to a dual optimization problem. While this approach may reduce the risk of negative estimates at $q_i$, it does not guarantee non-negativity of intensity elsewhere and undermines the computational advantages inherent in K$^2$IE due to the added complexity of dual optimization.

Finally, we discuss whether the problem in Eq. (11) is truly improper/wrong. (Kim, NeurIPS2024) showed that when RKHS kernels belong to the class of inverse M-kernels (IMKs), the corresponding equivalent kernels $h(x,x')$ are non-negative. This suggests that Eq. (11) may not be inherently improper because K$^2$IE, a sum of equivalent kernels, is also non-negative under IMKs. In one-dimensional cases, the Laplace kernel is known to be an IMK, but no construction is known for IMKs in higher dimensions, posing an interesting challenge. We will include the discussion in the text.

Can the method be used to reconstruction of intensity at unobserved regions

Yes, K$^2$IE can reconstruct intensity at unobserved regions. At submission, we had assumed the intensity estimation at observed regions. However, in light of your insightful comment, we revisited the model and confirmed that K$^2$IE in Eq. (12) could accept inputs from unobserved regions as well. Accordingly, a minor revision to Theorem 1 is warranted to reflect this more general setting. Specifically, all integral operators and the squared RKHS norm should be defined over the full domain ($X$ -> $\mathbb{R}^d$), and Eq. (13) should be updated as $q^*(x,s) = \delta(x-s) 1(s \in X) +\frac{1}{\gamma}k^*(x,s)$, where $1(\cdot)$ is the indicator. With this modification, Theorem 1, initially defined in $x \in X$, should be revised to hold for $x \in \mathbb{R}^d$. We sincerely appreciate the valuable comment. Due to time constraints, we have not yet conducted experiments to evaluate accuracy in unobserved regions. But, we are willing to include the results in the final version of the paper if you consider it essential.

How does the method compare to log-gaussian-cox process methods?

In response to several reviewers' suggestions, we conducted an additional experiment using the 2D synthetic dataset $\lambda_{\text{2D}}^{1.0}$ to include the result of a scalable Bayesian method​. This time, we adopted a variational Bayesian model with a quadratic link function (Lloyd, ICML2015), where 10 x 10 inducing points were employed. We appreciate your kind suggestion of references (we will cite them), but we have not yet become proficient with the R package. $L^2$, $|L|$, and cpu of the Bayesian method were 58.0(8.13), 5.61(0.47), and 28.4(0.96), respectively, where standard deviations are in brackets. The result highlights the high efficiency of K$^2$IE. Note that the reported results are based on 10 trials. We will include the results based on 100 trials in the text.

Does the methods provide uncertainty estimates for the reconstruction

K$^2$IE does not provide uncertainty estimates, which limits our model compared to Bayesian models.